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ABSTRACT 

This  thesis  concerns  itself  with  the  optimum-control  study 
of  a  class  of  infectious  diseases  which  are  propagated  through  the 
transfer  of  disease  micro-organism,  by  direct  contact  between  diseased 
and  healthy  individuals.  An  improved  model,  generally  applicable  to 
most  such  diseases,  reasonably  accurate  in  representing  the  phenomena 
connected  with  their  spread,  mathematically  simple,  and  suitable 
from  control  point  of  view,  has  been  formulated.  The  new  model,  while 
remaining  essentially  deterministic  in  nature,  incorporates  a  very 
desirable  feature  of  stochastic  models  by  representing  the  latent 
and  infectious  periods  of  disease  by  their  mean  values  and  respective 
standard  deviations. 

The  resulting  model  consists  of  a  set  of  non-linear 
differential  equations  which  include  functions  of  present  values  and 
past  history  of  both  state  and  control  variables.  Numerical  solutions 
have  been  obtained  for  the  simulated  model,  consisting  of  equivalent 
differential  difference  equations,  with  different  sets  of  parameters, 
but  without  control;  thus  demonstrating  the  applicability  of  the  model 
to  various  diseases.  The  effect  of  application  and  variation  of 
each  of  the  active  and  passive  vaccination  controls  has  also  been  studied. 

The  optimum  control  theory  has  been  applied  to  the 
problem  of  finding  the  most  economic  use  of  active  and  passive 


immunization  controls.  Application  of  Pontryagin's  Minimum  Principle 
to  this  case,  involving  functions  of  both  delayed  state  and  delayed 
control  ,  has  been  demonstrated  and  a  procedure  has  been  developed  for 
the  numerical  solution  of  the  resulting  control  problem. 

Using  the  numerical  procedure,  optimum  control  strategies 
have  been  obtained  for  different  values  of  reported  case  cost,  a 
parameter  representing  the  personal ,  social  and  economic  damage  done 
by  one  active  case  of  the  disease.  The  influence  of  delay  in  the 
effectiveness  of  active  control  on  the  resulting  optimum  cost  and 
controls  has  also  been  studied. 
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CHAPTER  I 

GENERAL  INTRODUCTION 

Z.  Z  Introduction 

Research  conducted  during  the  last  sixty  or  so  years 

in  the  fields  of  clinical  medicine,  biology,  epidemiology,  and  other 

related  disciplines  has  provided  considerable  insight  into  the 

phenomena  connected  with  the  spread  of  infectious  diseases.  A  large 

number  of  text  books  on  public  health,  containing  this  information, 

may  be  classed  as  non-technical  from  the  point  of  view  of  medical 

★ 

sciences;  references  [11],  [23],  [26],  and  [30]  are  a  few  examples  of 
this  fact.  This  easy  access  to  information,  coupled  with  the  common 
concern  of  every  individual  for  the  scourge  of  epidemics  to  which 
humanity  has  been  subjected  from  time  immemorial ,  has  attracted  the 
attention  of  a  considerable  number  of  researchers  from  sciences 
other  than  medical. 

Mechanisms  of  contagion  can  be  readily  classified  into 
two  sub-sections,  namely,  micro-mechanisms  and  macro-mechanisms;  the 
former  dealing  with  characteristic  properties  of  pathogenic  material 
of  infection  and  its  effects  on  the  biological  system  of  the  victim, 
and  the  latter  concerning  itself  with  the  behaviour  of  the  intended 
victims  as  a  group  and  the  susceptibilities  and  defenses  of  this 

*  Numbers  in  rectangular  brackets  indicate  references  listed  at  the 
end  of  this  thesis. 
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group  as  a  crowd  or  herd.  The  former  is  the  area  of  major  concern  of 
only  the  medical  or  biological  scientists,  whereas  the  latter  can  and 
has  benefitted  most  from  the  specialised  knowledge  of  mathematicians, 
statisticians,  ecologists,  and  system  scientists.  The  research 
reported  in  this  dissertation  falls  in  the  second  category. 

A  beginning  in  the  application  of  mathematical  concepts 
to  epidemics  may  be  said  to  have  been  made  with  the  studies  of  London 
Bills  of  Mortality  by  John  Graunt  and  William  Petty ,in  the  middle  of 
the  17th  century.  It  was  then  that  the  term  medical  statistics  really 
originated.  Early  attempts  on  the  study  of  statistical  returns  of 
deaths  from  smallpox  were  aimed  at  finding  some  empirical  laws  under¬ 
lying  the  spread  of  epidemics.  William  Farr,  in  fact,  succeeded  in 
fitting  normal  curves  to  quarterly  data  on  smallpox  deaths,  in  the 
middle  of  the  nineteenth  century.  Successes  like  this,  coupled  with 
research  in  the  biological  field  around  the  same  time,  laid  the 
foundations  of  the  present  mathematical  theory  of  epidemics. 

A  large  number  of  workers  have  proposed  a  number  of  math¬ 
ematical  models  to  mimic  the  behaviour  of  epidemics.  The  work  done 
in  this  field  has  been  extensively  reviewed,  and  good  number  of  review 
articles  and  publications  are  available;  prominent  among  them  are  those 
of  Serfling  [37],  Bailey  [2],  and  Deitz  [10],  Bailey's  monograph  [2] 
reviews  in  detail  the  work  done  up  to  1957,  and  has  been  updated  by 
Deitz  [10].  The  major  milestones  in  the  development  of  this  theory  are 
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the  works  of  Hamer,  Brownlee,  Lotka,  Ross,  Reed  and  Frost,  Kermick 
and  McKendrick,  and  E.B.  Wilson.  While  a  more  detailed  historical 
review  of  the  subject  is  postponed  to  the  later  parts  of  the  thesis, 
it  will  be  useful  to  discuss  here,  in  brief,  the  basic  concepts  of 
the  epidemic  theory  as  they  are  understood  today. 

1.2  Basic  Concepts 

Micro-organism  of  an  infectious  disease  is  carried  from 
a  diseased  individual  or  animal  to  the  healthy  one  either  by  direct 
contact  and  proximity,  or  through  an  intermediate  agent.  The  inter¬ 
mediate  agent  may  be  an  item  of  food,  water,  or  a  rodent  or  arthropod 
vector.  Attention  is  here  particularly  restricted  to  those  communicable 
diseases  which  pass  on  by  direct  contact  or  proximity.  The  causative 
agents  of  these  diseases  are  bacteria  or  virus.  It  may  be  pointed  out 
here,  however,  that  the  methods  discussed  in  this  thesis  may  be 
extended  to  other  classes  of  infectious  diseases  by  suitable  modi¬ 
fications. 

For  the  purpose  of  epidemic  study  a  population,  at  any 
time,  may  be  considered  to  consist  of  the  following  three  sub-populations 
susceptibleSj  infectives  and  immune s .  A  brief  description  of  these 
three  terms  is  in  order  here.  Swceptibles,  as  the  name  implies, 
are  the  healthy  individuals  who  have  a  very  weak  or  no  resistance  to 
the  disease  in  question.  Since  in  most  cases  an  individual  acquires 
immunity  to  a  disease  after  recovery,  it  may  be  assumed  that  a  member 
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of  the  susceptible  population  is  one  who  has  either  no  experience  of 
the  disease  in  question  or  has  lost  his  immunity  obtained  from  earlier 
exposures.  Infectives  are  the  individuals  who  have  sufficiently 
developed  pathogenic  material  in  their  system  so  that  they  are  capable 
of  passing  it  on  to  their  contacts,  thus  infecting  the  susceptible 
ones.  Immunes  are  the  members  of  population  who  have  acquired  immunity 
to  the  disease  in  question,  either  by  earlier  exposure  or  by  artificial 
means  of  immunization.  These  members  do  not  take  any  part  in  the 
process  of  disease  spread.  Therefore  the  individuals,  once  affected 
by  the  disease  but  since  dead,  recovered,  or  removed  from  circulation; 
either  by  hospitalization  or  effective  quarantine  on  the  appearance  of 
symptoms,  can  also  be  considered  in  the  immune  class. 

When  an  infective  somehow  gets  introduced  to  a  population 
of  susceptibles  and  immunes  only,  the  contacts  of  the  infective  pick 
up  the  disease  micro-organism  which  then  meets  the  biological  and 
chemical  defenses  of  the  host.  Depending  upon  the  degree  of  this 
resistance  and  its  own  capability  for  multiplication,  the  parasite 
establishes  itself  in  the  system  of  a  susceptible  host  after  a  period 
called  latent  period.  The  susceptible  now  himself  becomes  an  in¬ 
fective  and  starts  spreading  the  disease  to  his  healthy  contacts  till 
he  shows  the  symptoms  of  disease  and  is  removed  from  circulation  by 
quarantine,  hospitalization,  death,  or  recovery.  If  recovered,  he 
acquires  immunity  to  the  disease,  at  least  for  the  immediate  future, 
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and  enters  the  immune  population.  Some  relapses  in  diseases  like 
tuberculosis  do,  however,  occur.  Also,  some  infectives  develop  a 
tolerance  to  the  parasite  and  become  capable  of  passing  the  disease, 
at  least  for  a  limited  time,  without  themselves  ever  showing  the 
symptoms  of  disease;  they  are  called  carriers.  There  are  others 
who  become  immune  without  ever  becoming  sick  by  virtue  of  receiving 
repeated  but  small, doses  of  infection. 

The  time  elapsed  between  the  receipt  of  infection  by  a 
susceptible  and  his  becoming  infective  is  called  the  latent  period , 
and  that  between  first  exposure  and  appearance  of  symptoms  is  called 
the  incubation  period.  These  two  periods  have  definite  mean  values 
for  a  particular  disease.  The  difference  between  the  incubation 
period  and  latent  period  is  the  time  during  which  a  case  is  actively 
infective  and  is  hence  called  the  infectious  period. 

1.  3  Modeling 

The  discussion  so  far  leads  to  the  conclusion  that  the 
dynamics  of  a  crowd  disease  or  an  epidemic  can  be  represented,  at 
any  time,  by  the  numbers  of  susceptibles ,  infectives  and  recovered  or 
removed  individuals  and  the  rates  at  which  they  move  from  one 
category  to  the  other.  In  other  words,  an  epidemic  can  be  represented 
by  a  mathematical  model,  provided  parameters  for  the  rates  of  change 
of  the  above  variables  are  known.  In  fact,  a  number  of  such  models 
have  been  formulated  during  the  last  60  years,  as  shall  be  discussed 
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in  the  next  chapter.  The  primary  parameters  mentioned  above  are  the 
rates  at  which  effective  contacts  are  made  and  infectives  are  removed 
from  circulation.  However,  for  a  model  to  be  more  accurate,  in  addition 
to  the  consideration  of  chances  of  contact  and  removal,  the  mean  values 
and  the  upper  and  lower  limits  of  the  latent  and  incubation  periods 
should  also  be  incorporated.  Thus,  mathematically,  an  epidemic  can 
be  specified  by  the  effective  rate  of  contact,  and  mean  values  and 
variance  of  the  latent  and  incubation  periods. 

Since  the  events  of  contact  between  a  susceptible  and  an 
infective  (as  well  as  various  other  phenomena  connected  with  the  disease 
spread)  are  probabilistic  in  nature,  it  is  clear  that  the  model  re¬ 
presenting  an  epidemic  should  ideally  be  stochastic  in  nature.  However, 
when  dealing  with  large  populations,  a  deterministic  model  which  by 
its  nature  tends  to  evaluate  the  mean  values  of  variable  at  any  time 
rather  than  the  probabilities  of  these  numbers  is  sufficiently  accurate 
and  mathematically  much  simpler  than  its  stochastic  counterpart.  In  fact, 
a  stochastic  model  for  a  small  community,  when  considered  for  households 
of  more  than  two,  is  almost  unmanageable  mathematically.  Therefore, 
when  one  is  not  dealing  with  very  small  groups,  the  usefulness  of 
deterministic  models  can  not  be  over  stated. 

As  pointed  out  earlier,  models  help  in  understanding  the 
mechanism  of  the  spread  of  epidemics  at  macro  level.  The  understanding 
thus  gained  has  always  been  used  for  a  better  control  of  the  epidemics. 
Most  of  the  models  discussed  in  literature,  however,  can  be  said  to 
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be  predictive  in  nature  rather  than  control  oriented,  since  they  do 
not  include  parameters  representing  control  actions.  Because  of 
the  availability  of  methods  for  creating  artificial  immunity  to 
disease  by  use  of  active  and  passive  immunization,  and  because  of  the 
proven  effectiveness  of  chemoprophylaxis  in  preventing  disease, 
the  important  question  of  selecting  the  best  means  of  control  naturally 
arises.  If  only  the  predictive  models  without  control  parameters 
are  used,  the  control  strategy  arrived  at  can,  at  best,  be  said  to 
be  a  qualified  guess.  Whereas,  with  the  development  of  control  models 
and  application  of  optimization  techniques, a  more  logical  control 
strategy  can  be  evolved.  The  question  of  optimum  control  has  recently 
started  appearing  in  the  literature,  as  reviewed  in  subsequent  chapters. 

In  fact,other  attempts  on  some  kind  of  optimization  of  epidemic 
control  have  been  reported.  However,  to  the  best  of  this  author's 
knowledge  no  model,  sufficiently  accurate  and  yet  general  enough 
for  application  to  the  control  of  a  class  of  diseases,  has  been  reported. 
Nor  has  a  thorough  application  of  the  M inimum  P rinciple  of  optimum 
control  theory  to  this  class  of  problem  appeared  in  the  literature. 

An  improved  model  for  the  spread  of  contagious  diseases 
is  presented  in  the  next  Chapter.  The  new  model  is  basically  deterministic 
in  formulation  and  is  derived  from  the  standard  Kermack  and  McKendrick 
model.  However,  latent  and  incubation  periods  and  their  variation  for 
the  disease  in  question  have  been  incorporated  in  the  model.  The 
model  assumes  each  of  these  periods  to  be  normally  distributed  with 
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known  mean  value  and  standard  deviation,  and  these  values  constitute 
important  parameters  of  the  model.  Nominal  controls  in  the  form  of 
active  and  passive  immunizations  have  also  been  considered.  The 
resulting  model  has  been  used  to  evolve  a  most  economic  control 
strategy. 

Chapters  II,  III  and  IV  contain  the  major  contribution  of 
this  thesis.  Chapter  II  describes  the  new  model  and  its  parameters  in 
detail.  Numerical  solutions  of  the  model  equations,  presented  in 
this  Chapter,  clearly  show  its  advantages  over  the  earlier  models. 

The  Chapter  also  demonstrates  the  ease  with  which  the  model  can  be 
used  to  represent  various  diseases  by  adjusting  the  appropriate 
parameters.  The  introduction  of  control  and  its  effect  on  the  course 
of  disease  is  discussed  in  Chapter  III.  Two  different  controls, 
active  immunization  and  passive  immunization,  are  separately  considered, 
and  numerical  solution  for  various  predetermined  controls  are  pre¬ 
sented.  Chapter  IV  deals  with  the  extension  of  optimal  control 
theory  for  the  determination  of  best  control  strategy  between  the  two 
controls  used  simultaneously.  The  procedure  developed  minimizes  a 
cost  function  consisting  of  the  sum  of  the  cost  of  each  control  and 
that  of  each  case  affected  by  disease. 

It  is  hoped  that  this  dissertation,  in  addition  to  pre¬ 
senting  some  new  results  in  the  field  of  optimum  control  of  epidemics, 
also  represents  an  extension  and  good  example  of  some  recent  techniques 
in  optimum  control  theory.  The  new  model,  incorporating  the  two 
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different  controls  as  arrived  at  in  Chapter  III,  has  a  finite  number 
of  delays  in  both  state  and  control  variables.  Although  considerable 
work  on  optimum  control  of  systems  with  delays  in  state  or  control  has 
been  reported,  and  recently  some  papers  discussing  delays  both  in 
state  and  control  have  also  appeared,  yet  the  author  has  not  come 
across  any  example  of  such  an  application  in  the  literature.  Numerical 
solutions  presented  in  this  thesis  will,  hopefully,  provide  a  good 
example  of  practical  application  of  optimum  control  to  a  nonlinear 
system  with  constraints,  and  having  finite  number  of  delays  in  both 
state  and  control . 
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CHAPTER  II 

NEW  MODEL 

2. 1  Introduction 

Introductory  remarks  made  in  the  last  chapter  suggest 
the  existance  of  a  large  number  of  mathematical  models  suitable  for 
representing  the  epidemic  spread.  In  fact,  the  very  understanding 
of  the  mechanisms  of  disease  spread,  which  we  now  possess,  is  mainly 
due  to  the  liberal  use  of  these  models  during  the  last  half  century. 

Among  the  available  models,  both  deterministic  and  stochastic,  some 
are  specific  to  special  cases  whereas  the  others  are  more  general. 

Again,  some  are  relatively  simple  whereas  the  others  are  more  complicated. 
Generally  speaking,  simple  ones  are  less  accurate  in  the  representation 
of  epidemic  spread.  Therefore  a  model  which  is  accurate  enough  and 
yet  mathematically  simple  is  always  desirable.  One  such  model  is 
introduced  in  this  chapter. 

The  new  model  is  an  improved  version  of  a  well  known 
deterministic  model  and  incorporates  many  desirable  features  of  other 
models,  both  deterministic  and  stochastic.  This  model  is  quite  general 
in  application  to  various  diseases  and  quite  suitable  for  introducing 
control  parameters.  Moreover,  it  is  also  suitable  for  the  optimum 
control  analysis. 


-  11 


The  new  model  is  discussed  in  detail  in  this  chapter, 
whereas  the  analysis  of  the  model  with  respect  to  controls  is  taken  up 
in  the  subsequent  chapters. 

2.  2  Historical  Review 

The  extensive  literature  on  the  mathematical  theory  of 
epidemics  has  been  adequately  reviewed  from  time  to  time.  Serfling  [37] 
presented  a  comprehensive  review  of  the  state  of  the  art  as  it  existed 
in  1952.  In  1957,  Bailey  [2]  published  his  now  famous  monograph 
summarising  all  the  work  in  the  field  up  to  that  time.  This  work  did 
more  for  the  field  of  mathematical  modeling  of  epidemics  than  any  other 
single  effort.  It  is  no  surprise,  therefore,  that  Bailey  predicted  the 
future  use  of  mathematical  models  in  the  economical  control  of  epidemics, 
which  is  now  receiving  some  attention  in  the  literature  and  is  also 
the  subject  of  this  author's  present  effort.  Bailey's  work  has  been 
adequately  supplemented  by  the  review  article  of  Deitz  [10],  which  up¬ 
dates  it  to  1967.  These  three  efforts  alone,  not  withstanding  the 
parallel  efforts  of  other  authors,  present  an  almost  complete  picture 
of  the  state  of  the  art  in  the  field.  Any  attempt  to  improve  upon  this 
extensive  coverage  of  historical  reviews  is  neither  possible  nor  desirable 
here.  Some  general  comments,  necessary  for  the  understanding  of  the 
new  model,  are,  however,  made  in  the  next  few  paragraphs. 

The  possibility  of  defining  some  empirical  relations  describ¬ 
ing  the  epidemic  spread  was  first  visualized  by  William  Farr  in  1840. 
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As  pointed  out  in  the  last  chapter,  early  efforts,  including  those  of 
Farr,  were  devoted  to  the  finding  out  of  these  relations  from  the 
published  statistical  case-report  records  by  the  use  of  mathematical 
curve  fitting  methods.  The  works  of  Hamer,  Brownlee,  Lotka  and  other 
pioneers  in  the  field,  were  also  directed  to  the  same  goal.  These 
investigations  helped  establish  the  fact  that  the  number  of  cases 
reported,  at  any  time,  depends  upon  the  number  of  susceptibles  and 
infectives  in  the  population.  This  line  was  further  advanced  by  Ross, 
Reed  and  Frost,  Wilson  et  al  and  the  Kermack-McKen'drick  team,  among 
others.  Most  of  the  early  attempts  were  deterministic  in  nature  for  the 
sake  of  mathematical  simplicity.  The  most  popular  deterministic  model 
is  the  one  by  Kermack  and  McKendrick  [21],  proposed  in  1927.  It  was 
on  the  basis  of  this  model  that  the  authors  proposed  the  now  celebrated 
Threshold  Theorem  which  states  that  for  a  disease  to  start  the  pop¬ 
ulation  of  susceptibles  should  be  more  than  a  certain  threshold  value. 
This  model  is  the  one  most  often  used,  and  it  forms  the  basis  of  the 
model  to  be  proposed  in  this  chapter. 

Dissatisfaction  with  deterministic  models,  due  to  their 
inaccuracy  in  representing  essentially  probabilistic  phenomena,  led  to 
the  development  of  mathematically  more  complicated  stochastic  models. 

The  simple  stochastic  epidemic  model,  which  led  to  those  used  today, 
is  due  to  Bartlett  [6]  and  Bailey  [2],  [3],  [4],  [5],  and  was  first 
introduced  by  Bartlett  in  1949.  The  accurate  consideration  of  latent 
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and  incubation  periods  was  first  introduced  in  the  stochastic  models, 
and  some  concepts  evolved  there  are  used  in  the  model  proposed  in  this 
chapter.  Another  significant  but  largely  unexplored  area  is  the 
geographical  element  introduced  in  epidemic  models,  mainly  by  Menther 
and  Bailey  [2] ,  [4] . 

To  sum  up,  two  kinds  of  epidemic  models  are  now  available 
in  the  literature  :  deterministic  and  stochastic.  The  standard 
deterministic  model  is  the  Kermack  and  McKendric  [21]  model.  The 
stochastic  models,  however,  are  being  continuously  improved  upon  as 
increasing  numbers  of  researchers  adopt  them  because  of  their  better 
accuracy  and  better  representation  of  the  phenomena  connected  with 
disease  spread.  The  immediate  appeal  of  the  stochastic  models  stems 
from  the  fact  that  everything  connected  with  disease  spread  depends 
upon  chance,  and  these  models  take  this  fact  into  account  precisely. 

But  when  large  populations  are  involved,  it  follows  from  the  low  of 
large  numbers  that  the  stochastic  deviations  are  small  compared  to  the 
average  values,  hence  the  argument  that  deterministic  models  are  useful 
can  not  be  rejected  out  of  hand.  This  is  especially  so  when  the 
mathematical  simplicity  of  these  models  far  outweighs  the  marginal 
advantages  of  stochastic  models  when  applied  to  larger  populations.  In 
fact,  due  to  their  mathematical  simplicity,  the  deterministic  models 
are  the  only  ones  suitable  for  effective  optimum  control  studies  at 
the  present  state  of  development. 
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It  is  due  to  the  above  considerations  that  the  model  to 
be  introduced  is  kept  essentially  deterministic,  even  though  the  latent 
and  incubation  or  infectious  periods  and  their  statistical  variations 
are  accounted  for.  The  new  model  is  quite  general  in  application  and 
suitable  from  the  control  point  of  view.  Since  the  new  model  is  based 
on  the  existing  Kermack  and  McKendrick  model  and  is  an  improvement  on 
it,  it  is  proposed  to  discuss  the  old  model,  in  some  detail,  in  the 
next  section. 

2.3  Kermack  and  McKendrick  Model 

The  standard  deterministic  model  for  a  closed  population 
proposed  by  Kermack  and  McKendrick  and  mentioned  in  the  last  section  is: 


dx 

dt 

My 

(2.1) 

d Z  = 

3  x  y  -  y  y 

(2.2) 

dt 

dz  = 

Y  y 

(2.3) 

dt 

where  x,  y,  z  are, 

respectively,  the  numbers  of  susceptibles. 

infecti ves 

and  removed  or  recovered  cases  at  any  time,  3  is  the  contact  rate  for 
the  disease  in  question  and  y  is  the  removal  rate.  The  solution  of  this 
model  gives  a  bell  shaped  case  report  curve  of  the  type  generally 
observed  in  actual  practice.  This  model  is  a  landmark  in  the  history  of 
the  mathematical  theory  of  epidemics,  since  the  now  celebrated  threshold 
theorem,  that  an  epidemic  will  not  start  until  and  unless  the  number 
of  susceptibles  is  more  than  the  relative  removal  rate  p=y/3,  was 
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derived  from  this  model.  In  view  of  the  non-1 inearity  in  the  differential 
equations  epidemiological  conclusions,  like  the  threshold  theorem,  were 
drawn  from  equilibrium  conditions  [29],  and  results  were  compared  with 
those  of  a  disease  in  rat  populations  under  laboratory  conditions. 

It  may  be  pointed  out  here  that  the  later  stochastic  models  also  arrived 
at  a  threshold  theorem  similar  to  that  mentioned  above,  thus  proving 
the  essential  correctness  of  the  model. 

Although  mathematically  very  simple,  the  model  is  not  very 
accurate  in  representing  epidemics,  especially  in  the  beginning  and  near 
the  end  of  the  disease  cycle.  This  model  assumes  the  removals  (yy)  to 
be  proportional  to  the  population  of  infectives.  This  assumption  is 
very  inaccurate,  especially  at  the  terminals  of  the  curve.  The  model, 
by  virtue  of  its  formulation,  assumes  a  zero  latent  period  and  a  negatively 
distributed  (Poisson)  infectious  period  [21];  this  may  be  somewhat 
correct  for  quick  spreading  diseases  like  influenza  and  common  cold, 
but  not  so  for  most  other  diseases. 

Due  to  its  mathematical  simplicity  this  model  has  been  used 
extensively,  and  it  is  also  the  basis  of  ReVelle's  model  for  economical 
control  of  T.B.  [33],  [34],  [23].  His  model  consists  of  nine  differential 
equations  instead  of  the  three  of  the  Kermack  McKendrick  model.  The 
formulation  of  these  equations  in  the  two  models  is  very  similar.  The 
various  subpopulations  considered  in  ReVelle's  tuberculosis  model  are 
those  of  susceptibles ,  susceptibles  vaccinated  with  B.C.G.,  infected 
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but  non-active  cases,  non-actives  with  prophylaxis  cured  cases,  and 
naturally  recovered  cases.  This  improved  model  appears  to  be  re¬ 
asonably  accurate  for  application  to  the  control  of  tuberculosis,  but 
still  does  not  consider  latent  and  infectious  periods  explicitly.  It 
also  uses  the  proportional  rates  of  the  type  'y 1  used  in  equation  (2.2) 
to  represent  transfers  of  population  from  one  group  to  the  other  during 
disease  spread.  Since  this  model  is  highly  specialized  for  tuberculosis, 
its  application  to  other  situations  would  be  very  complicated,  if  not 
impossible. 

2.  4  New  Model 

Reference  to  standard  text  books  on  public  health  :  Leavell 
et  al  [26],  and  epidemiology  :  Taylor  et  al [40] ,  shows  that  the  in¬ 
cubation  period  for  every  disease  varies  between  a  lower  and  an  upper 
limit  and  these  values  of  incubation  period  are  so  consistent  that 
they  are  used  for  differential  diagnosis  of  diseases.  This  fact  is 
further  borne  out  by  the  published  results  of  studies  conducted  by  Hope 
Simpson,  [14],  [15].  Reference  [15]  tabulates  the  distribution  of  cases 
against  incubation  period  in  days  for  three  different  diseases,  namely  : 
measles,  mumps  and  chickenpox  (varicella).  Bailey  has  discussed  the 
Hope  Simpson  data,  in  some  detail,  in  chapter  seven  of  reference  [2], 
and  has  concluded  that  normal  distribution  curves  can  be  fitted  to  the 
incubation  period  data.  Continuing  his  discussion,  on  the  basis  of  his 
earlier  published  work,  he  further  concludes  that  it  is  reasonable  to 
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assume  a  normally  distributed  latent  period  and  a  constant  infectious 
period.  The  fact  that  incubation  period  and  latent  period  have  some 
kind  of  distribution  was  also  proved  by  Abbey  [1]  and  Sartwell  [35]. 

Although  some  earlier  deterministic  models  [44]  used 
constant  latent  and  infectious  periods,  yet  it  is  only  the  recent 
stochastic  models  which  use  the  variations  in  these  periods  as  well. 
Bailey  [2]  proposed  the  stochastic  model  considering  normally  distributed 
latent  period  and  constant  infectious  period  (which  implies  a  normally 
distributed  incubation  period,  since  incubation  period  is  the  sum  of 
latent  and  infectious  periods).  He  outlined  a  procedure  for  the 
estimation  of  parameters  for  this  model  and  calculated  these  parameters 
for  an  epidemic  of  measles  using  Hope  Simpson  [15]  data.  Among  the 
parameters  calculated  were  the  mean  value  and  standard  deviation  of 
the  latent  period  and  the  value  of  infectious  period  (assumed  constant 
in  his  model).  Significantly,  Bailey  and  Stenberger  [5]  have  recently 
(1970)  revised  the  above  procedure  for  estimation  of  parameters  for  the 
sake  of  better  accuracy.  In  addition  to  the  recalculation  of  parameters 
for  measles  they  have  also  calculated  parameters  for  infectious  hepatitis 
using  the  recent  (1959-67)  data  privately  supplied  by  Dr.  K.  Peterson  of 
Hamburg. 

It  is  clear  from  the  above  discussion  that  incubation  period 
for  every  disease  has  a  distinct  mean  value  and  a  definite  distribution. 
Although  incubation  period  is  directly  observable  from  the  case  histories 


-  18  - 


of  the  patients,  the  latent  period  can  only  be  found  indirectly. 

Bailey's  assumption  of  a  constant  infectious  period,  while  considering 
the  latent  period  to  be  normally  distributed,  seems  to  have  been 
made  at  least  as  much  for  mathematical  simplicity  of  his  stochastic 
model  as  for  accuracy.  In  the  opinion  of  this  author,  a  more  general 
model  should  have  all  the  periods  represented  by  their  mean  values  and 
standard  deviations.  Consideration  of  a  distributed  rather  than  a 
constant  infectious  period  may  also  be  made  to  compensate  for  the  fact 
that  many  reported  cases  may  escape  effective  removal  from  circulation, 
and  may  be  spreading  the  disease,  at  least  for  some  time,  after  the 
appearance  of  symptoms. 

We  can  now  proceed  to  incorporate  the  above  feature  in  the 
epidemic  model.  It  may  be  pointed  out  here  that  for  an  accurate  model 
we  need  the  mean  value  and  standard  deviations  of  latent  period  and 
those  for  one  of  the  other  two  periods,  i.e.  either  infectious  period  or 
the  incubation  period.  A  version  of  the  new  model  may  be  stated  as 
follows  : 


X1  ■ 

-  eo  X1  X2 

(2.4) 

II 

CXI 
•  X 

+-> 

CC 

1 

+-> 

<c 

(2.5) 

II 

CO 
•  X 

R(t) 

(2.6) 

II 

•  X 

K(t)  B0  X1  X2  , 

(2.7) 

and 

X,  are  the  numbers  of  susceptibles. 

infectives  and 
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removed  or  recovered  cases,  respectively,  and  ^ ,  X2,  X3  are  the 
corresponding  rates  of  change.  is  the  rate  at  which  infected  (but 
not  yet  infective)  cases  are  generated;  $Q  is  the  effective  contact 
rate;  K(t)  is  the  fraction  of  population  which  finally  becomes  active 
in  spreading  the  disease,  [1  - K( t ) ]  being  that  fraction  which  becomes 
immune  due  to  repeated  but  small  doses  of  infection  [40,  page  142]. 

A(t)  is  the  rate  at  which  infected  cases  become  infectives  after  their 
latent  period  is  over,  and  R(t)  is  the  rate  at  which  infectives  show 
symptoms  and  are  reported,  and  hence  are  considered  removed  from 
circulation.  Expression  for  variables  A(t)  and  R(t)  are  derived  in 
the  next  few  paragraphs  and  are  represented  by  equations  (2.8)  to  (2.14). 

The  model  is  clearly  based  on  the  earlier  Kermack  McKendrick 
model.  However  the  manner  of  representing  the  introduction  of  active 
cases  and  removal  of  reported  cases  is  different  in  the  two  models. 

It  is  now  possible  to  calculate  the  active  cases  and  reported  cases  on 
the  basis  of  latent  and  incubation  periods.  The  expressions  for  these 
calculations  are  discussed  next. 

It  was  concluded  in  the  earlier  discussion  that  it  is 
reasonably  accurate  to  assume  the  latent,  infectious  and  incubation 
periods  to  be  normally  distributed.  This  hypothesis  will  now  be  used  to 
evaluate  the  variables  A(t),  and  R(t).  We  refer  to  figure  2.1  for 
deriving  the  expression  for  the  calculation  of  A(t).  Part  (a)  of  this 
figure  shows  a  probable  curve  for  the  past  history  of  infection  rate  X^ 
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FIGURE  2-1 


and  part  (b)  represents  a  normally  distributed  unit  impulse  with  its 
distribution  parameters  and  a-j  being  the  mean  value  and  standard 
deviation  of  the  latent  period.  The  two  curves  are  drawn  on  the  same 
time  axis.  The  number  of  cases  becoming  active  at  a  time  t,  from  among 
those  who  picked  up  infection  x  days  back,  will  be  the  product  of  the 
ordinates  of  two  figures  at  time  t-x.  Mathematically  this  number  will 
be 


1 

y  2tt 


exp 


1 


X^(t-x) , 


and  the  total  number  of  newly  active  infectives,  appearing  at  time  t  will 
be  the  integral  of  this  product  over  x  varying  from  0  to  °o.  Therefore 


A  ( t )  =  -L 

y  2tt 


X4(t-x) 


•  exp 


1 


(2.8) 


Since  the  ordinate  of  the  normal  distribution  curve  is  nearly  zero  for  T 
greater  than  2y^  equation  (2.8)  can  be  approximated  to 
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Similarly,  if  y^  and  are  the  corresponding  mean  value  and  standard 

deviation  of  incubation  period,  the  expression  for  the  reported  cases 

becomes , 
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R(t)  can  also  be  calculated  from  the  past  history  of  A(t)  rather  than 
that  of  X^.  The  correspond!' ng  expression  then  is. 


2y2 

R(t)  =  — —  —  f  A(t-x)  •  exp 

/  2tt  o0  J 

C  0 

if  u2  and  a 2  are  the  mean  values  and 
period  instead  of  being  those  of  the 


2a2 

standard  deviation  of 
incubation  period. 


(2.11) 

infecti ous 


Referring  again  to  figure  2.1,  we  can  rewrite  the  expressions 
for  A(t)  and  R(t)  in  discrete  form  by  making  the  valid  approximation  that 
at  any  time  t-nx,  the  area  under  the  curve  between  times  t-nx  and  t-(n-l)x, 
is  equal  to  the  ordinate  at  (t-nT)  multiplied  by  1,  where  n  is  a  number 
and  t  represents  an  interval  of  one  day.  This  area  is  shown  shaded  in 
the  figures  2.1  (a)  and  2.1  (b)  for  a  typical  period  of  one  day.  With 
this  assumption  the  expressions  for  A(t)  and  R(t),  in  discrete  form,  are: 
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when  incubation  period  is  used,  or 


R(t)  = 
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when  infectious  period  is  used. 
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Normalization 

Let  us  now  examine  the  model  for  dimensional  consistency. 

A  dimensional  analysis  of  equation  (2.4)  shows  that  for  X-j  ,  and  X-j 
to  represent  actual  numbers,  in  a  population,  bq  must  have  the 
dimension  of  inverse  population.  Thus  other  conditions  being  same,  the 
value  of  6q  will  be  different  for  communities  of  different  populations. 
The  following  example  illustrates  this  point  more  clearly. 

Suppose  a  community  of  10,000,  consisting  entirely  of 
suscepti bl es ,  has  10  infectives  introduced  at  a  given  time.  Using 
equation  (2.4)  and  a  contact  rate  e  -j ,  the  number  of  infected  cases 
generated  per  day  (since  day  is  the  unit  of  time  used  in  this  dissert¬ 
ation)  will  be 

Bol  x  10  x  10,000  =  105  &ol.  (2.15) 

If  we  now  consider  the  same  community  to  be  consisting  of  10  sub¬ 
communities,  each  of  1,000  susceptibl es ,  and  one  infective  is  introduced, 
simultaneously,  in  each  of  these  subcommunities,  the  total  number  of 
cases  generated  per  day  will  be 

10  (bo2  X  1  X  1,000)  =  104  eo2  (2.16) 

where  3  .  is  the  new  contact  rate  for  each  community.  Since  both  the 
o2 

above  expressions  represent  the  same  situation,  3  =  10  . 
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For  the  model  to  be  more  general,  however,  a  contact  rate 
independent  of  population  size  will  certainly  be  better.  This  was 
achieved  by  ReVelle  [”33],  in  his  T.B.  model,  by  representing  susceptibles 
as  a  fraction  of  total  population  but  maintaining  the  infectives  as 
actual  number;  the  infection  rate,  then,  was  also  in  actual  numbers.  A 
still  better  method  of  achieving  the  same  result,  in  the  opinion  of  this 
author,  is  to  normalize  all  the  state  (population)  variables,  i.e.  to 
represent  the  state  variables  as  fractions  of  total  population  rather 
than  actual  numbers. 

Therefore,  replacing  the  upper  case  variables  by  their  lower 
case  counterparts  so  that  the  new  variable  is  the  old  one  divided  by 
total  population  N,  i.e.  x  =  ^  ,  the  model  can  be  rewritten  as 
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-  3  x]x2 
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where  3  =  N3Q^ 


(2.21) 


-  25  - 
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if  infectious  period  is  used  and  reported  cases  are  calculated  from  the 
past  history  of  active  rate.  Corresponding  values  aft)  and  r(t)  in 
discrete  form  are 
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(2.25) 


(2.26) 


as  the  case  may  be. 

The  variables  x-j ,  x2»  and  x^  now  represent,  respectively,  the  sub¬ 
populations  of  susceptibles ,  infectives,  reported  cases  and  infected 
cases  as  fractions  of  the  total  population,  and  x-j ,  x^,  x^  and  x4 
represent  their  respective  rates  of  change,  aft)  and  r(t)  are  the 
corresponding  active  and  report  rates,  again  as  fractions  rather  than 
actual  numbers.  The  contact  rate  3  is  now  independent  of  the  population 


' 
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size,  thus  making  the  model  more  universal  in  application. 

The  term:  per  unit  value ,  commonly  used  in  electrical 
engineering  terminology,  can  be  applied,  successfully,  to  describe  the 
normalized  quantities,  defined  above.  Here,  as  there,  this  term  means 
that  a  quantity  called  per  unit  value  is  the  ratio  of  actual  value  to 
a  reference  value.  In  the  present  context,  if  x  =  ^  is  a  sub-population 
expressed  as  a  per  unit  value,  then  X  is  the  actual  sub-population  and 
N  is  the  total  population  or  reference  population.  Given  a  per  unit 
quantity,  an  actual  value  can  be  found  by  multiplying  per  unit  value  by 
the  corresponding  reference  value.  All  the  variables  represented  by 
lower  case  letters,  in  the  subsequent  discussion,  will  represent  per 
unit  quantities. 

Contact  Rate  3 

The  contact  rate  3,  used  in  the  normalized  new  model  has  the 
same  dimension  as  that  of  infection  rate  (3)  used  by  ReVelle  [33]. 
Therefore  both,  heuristic  and  probabilistic,  derivations  of  infection 
rates  given  by  ReVelle  [33,  pp.  26-30],  are  equally  applicable  in  the 
present  case  also.  Physical  interpretation  of  3  is  best  understood  by 
quoting  ReVelle  [33,  D.  26]. 

"3  =  average  number  of  individuals  per  unit  time  that  any 
individual  (active  case  or  not)  will  encounter 
sufficiently  to  cause  infection.  The  value  does  not 
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depend  on  whether  the  person  encountered  is  susceptible 
or  the  individual  doing  the  encountering  is  infectious. 
However,  only  if  the  encounterer  is  infectious  and  the 
encountered  susceptible  will  a  new  infection  arise.  An 
individual  encountered  sufficiently  is  call  a  "contact"; 
then  8  is  the  average  number  of  contacts  any  individual 
makes  per  unit  time.  The  parameter  8  is  characteristic 
of  the  disease  and  the  average  individual's  behaviour." 

Elaborating  the  last  sentence  of  the  above  quotation  we  can 
say  that  8  depends  upon  the  infectiousness  of  the  disease,  weather  and 
meteorological  conditions,  and  living,  working  and  social  conditions  of 
the  community.  A  practical  method  of  evaluating  8  is  by  curve  fitting 
on  available  data.  However,  with  sufficient  knowledge  of  the  quantitative 
effects  of  the  above  mentioned  social  and  biological  factors  and  with 
enough  experience,  an  a-priori  evaluation  of  8  is  possible,  at  least 
theoretical ly. 

2.5  Comments  on  the  model 

As  pointed  out  earlier,  the  important  features  of  the  new 
model  are:  (1)  normalization,  (2)  representation  of  the  latent,  infectious 
and  incubation  periods  by  their  mean  values  and  standard  deviations. 
Normalization  provides  a  better  interpretation  of  contact  rate  8  and 
makes  it  independent  of  the  size  of  the  population  to  which  the  model  is 
applied.  The  representation  of  sub-population  as  fraction  has  also  been 
used  to  advantage  by  Bailey  T3]  for  perturbation  approximation  to  simple 
stochastic  epidemics.  A  somewhat  similar  representation  has  also  been 
suggested  by  Landau  and  Rapoport  f24]  in  their  spread  of  information 
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model;  the  expression  arrived  at  there  is  somewhat  similar  to  the 
expression  we  would  obtain  after  substituting  a(t)  and  r(t)  in  equation 
(2.18).  The  time  ' t 1  in  our  model  is  similar  to  the  "private  time" 
used  in  that  paper. 

Hope  Simpson  [15],  in  addition  to  arriving  at  the  distrib¬ 
ution  of  the  incubation  period  for  three  different  diseases  as  discussed 
earlier,  also  showed  that  the  infectiousness  of  each  disease  is 
different.  Thus,  other  conditions  being  equal,  6  for  one  disease  can 
be  converted  to  that  for  the  other  by  multiplying  it  by  the  ratio  of 
infectiousness  of  the  two  diseases.  As  an  example,  Hope  Simpson  [15] 
calculated  the  infectiousness  of  measles,  varicella  and  mumps  to  be 
66.5%,  48.2%  and  32.1%.  Thus  if  6  for  mumps  is  1, 

6  for  varicella  =  1  x  ~  1.5 

32.1 

rr  r 

and  3  for  measles  =  1  x  — —  -  2.1 

32.1 

Thus  we  find  that,  as  far  as  the  model  is  concerned,  a 
disease  is  identified  by  3,  y-j  ,  y^,  a-j  and  c^,  and  these  constants  are 
unique  for  each  disease.  So  our  model  can  be  made  to  represent  any 
disease  for  which  these  constants  can  be  identified. 

The  new  model  has  available,  at  any  time,  the  values  of 
variables  representing  suscepti bl es ,  infectives,  infected  cases,  reported 
cases,  infected  case  rate  and  reported  case  rate.  Since  control  of  an 
epidemic  is  affected  by  modification  of  one  or  more  of  these  variables. 
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introduction  of  control  parameters  is  quite  easy.  Therefore  the  model 
is  quite  suitable  for  application  of  control.  This  aspect  will  be 
discussed  in  greater  detail  in  the  next  chapter. 

2.6  Analysis  of  the  model 

It  was  pointed  out  earlier  that  the  simple  model  of  Kermack 
and  McKendrick  already  presented  many  difficulties  for  its  mathematical 
solution.  The  modified  model,  as  is  apparent  from  equations  (2.17)  to 
(2.22),  is  much  more  complicated  because  of  the  integro-differential 
nature  of  the  equations.  Therefore  no  general,  closed-form  solution 
is  attempted.  The  model  is  analysed  here  by  digital  simulation  and  the 
important  results  are  presented. 


The  procedure  used  for  solving  the  model  numerically  is  to 
calculate  the  values  of  a(t)  and  r(t)  for  each  day  from  the  past  stored 
values  of  the  variables  representing  sub-populations  (now  called  state 
variables).  Expressions  (2.24)  to  (2.26)  are  used  for  calculating  the 
variables  numerically.  These  expressions  are  greatly  simplified  if  new 
weight  multiplier  vectors  are  defined  as  below: 


WLP(n)  =  — 

1_ 

exp 

1  ,  ,2 
-  2  (ur  nr) 

2ol  J 

/  2tt 

°1 

WIP(n)  =  — 

exp 

r  i  (  ,2  ’ 

L  z°2  J 

/  2t t 

°2 

(2.27) 

(2.28) 


Then 


a  ( t ) 


WLP(n)  •  x4(t-nr) 

n=o 


(2.29) 
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WIP(n)  •  x^(t-nx) 


(2.30) 


n=o 


for  incubation  period  or 


(2.31) 


n=o 

for  infectious  period. 

The  values  WLP  and  WIP  now  represent  the  ordinates  of  the 
normal  distribution  curve  at  each  day  between  the  current  time  and  2y 
days  back.  The  normal  distribution  is  nearly  zero  at  y  days  away  on 
either  side  of  the  mean  peak,  and 


^  WLP(n)  =  1 


WIP(n)  =  1. 


n=o 


Keeping  the  calculated  rates  a (t)  and  r(t)  constant  for 


one  day,  the  numerical  solution  of  the  differential  equations  is  obtained 
for  one  day  by  using  the  Runga  Kutta  procedure.  Thus  the  final  value  of 
variables  at  the  end  of  the  day  is  obtained  from  the  known  initial  values. 
The  final  values  now  calculated  constitute  the  initial  values  for  the 
next  day.  The  process  is  continued  till  the  value  of  infectives  becomes 
negl igibly  small . 
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The  solution  of  the  model  is  obtained,  taking  the  mean 
values  of  latent  and  infections  period  to  be  7  days.  The  standard 
deviations  assumed  are  1.4  and  2.0  respectively.  Although  these  values 
are  completely  arbitrary,  they  resemble  those  calculated  by  Bailey  [2] 
for  his  measles  epidemic  model.  (Since  the  purpose  here  is  only  to 
demonstrate  the  applicability  of  the  new  model  to  various  diseases,  the 
arbitrariness  of  the  above  values  does  not  detract  from  the  ensuing 
discussion.)  The  model  has  been  solved  for  values  of  e  varying  from 
0.2  to  3.0  and  for  initial  susceptible  populations  varying  from  0.1  to 
1.0.  Since  the  only  observable  measure  of  disease,  in  actual  practice, 
is  the  number  of  reported  cases  per  day  or  week,  the  results  of  the 
above  analysis  have  been  plotted  with  x-axis  representing  time  in  days 
and  y-axis  representing  the  reported  cases  (x^)  per  dav. 

Figure  2.2  shows  the  solution  of  the  new  model  with  a 
constant  3  but  with  variable  initial  suscepti bl es .  The  result  is  a 
family  of  10  bell  shaped  curves,  each  representing  the  response  for  one 
value  of  initial  susceptibles  between  0.1  to  1.0.  The  figure  shows 
clearly  that  the  spread  of  disease  is  most  violent  when  the  entire  popul¬ 
ation  is  susceptible,  i.e.  when  disease  has  been  introduced  into  a  virgin 
population.  The  severity  of  the  disease  successively  decreases  as  the 
initial  susceptibles  are  assumed  lesser  and  lesser,  till  a  threshold 
value  of  initial  susceptibles  is  reached  below  which  disease  does  not 
start.  The  latter  point  is  made  more  clear  in  figure  2.3,  which  is  a 
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replot,  at  enlarged  scale,  of  the  area  enclosed  bv  dotted  lines  on  figure 
2.2.  The  threshold  value  of  initial  susceptibles  is  0.1  for  6  =  1.5. 

These  results  are  in  line  with  the  known  results  of  standard  deterministic 
models. 

In  addition  to  the  above  conventional  results,  we  observe 
that  in  the  beginning  of  the  epidemic  the  successive  generations  of 
cases,  separated  by  approximately  an  incubation  period,  are  also  visible. 
This  fact,  generally  observed  during  actual  epidemics,  is  shown  very 
clearly  in  figure  2.3.  Therefore,  we  note  an  improvement  over  previous 
deterministic  models  which  failed  to  represent  this  phenomenon. 

Figure  2.4  shows  the  rate  of  reported  cases  as  function  of 
time,  now  with  variable  $  and  constant  initial  susceptibles.  6  is  varied 
over  a  wide  range  (0.2  to  3.0)  whereas  initial  susceptibles  are  the  same 
(0.5)  in  each  case.  We  find  that  for  a  given  initial  susceptibles  the 
epidemic  peak  is  larger  for  larger  values  of  3.  Figure  2.5,  which  is 
again  a  replot,  at  enlarged  scale,  of  the  area  enclosed  by  dotted  lines 
on  figure  2.4,  shows  that  there  is  a  threshold  value  of  3  below  which 
disease  will  not  start.  Thus  we  find  that  each  initial-susceptibles 
value  has  a  threshold  value  of  3.  This  result  is  in  line  with  the  mo¬ 
dified  threshold  theorem  of  Landau  et  al .  [24],  although  stated 
differently. 


The  existence  of  threshold  values  of  initial  susceptibles 


and  contact  rate  6  gives  a  possible  explanation  of  the  flare-up  of  endemic 
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FIG.  2-5  BOXED  AREA  IN  FIGURE  2-4 ON  ENLARGED  SCALE 
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diseases.  It  appears  that  when  an  epidemic  wanes  it  leaves  most  people 
immune  to  the  disease,  i.e.  the  "herd  immunity"  is  high.  As  time  passes 
the  herd  immunity  goes  on  decreasing,  due  to  the  addition  of  more 
susceptibles  by  new  births  and  removal  of  immune  older  people  by  death. 
Some  individuals  also  loose  immunity  with  time.  Thus  the  susceptible 
population  builds  up  but  still  remains  below  the  threshold  value  for 
the  prevailing  6.  But  when,  due  to  meteorological  changes  or  some 
other  reasons,  3  suddenly  drops,  the  initial  susceptibles  are  now  more 
than  the  threshold  value  for  the  new  6  and  a  flare-up  of  epidemic 
occurs.  Studies  of  Spicer  [39]  and  Wise  [46],  support  the  contention 
that  weather  changes  and  meteorological  factors  have  definite  effect  on 
the  infectiousness  of  diseases,  thus  supporting  the  above  hypothesis. 

In  conclusion,  we  can  say  that  the  new  model  is  quite 
reasonable  in  representing  many  phenomena  connected  with  disease  spread 
and  is  likely  to  be  a  useful  tool  in  the  study  of  epidemics.  The  latter 
is  borne  out  by  the  control  studies  to  be  presented  in  the  next  two 
chapters.  It  may  be  mentioned  here  that  the  curves  of  figure  2.4  are 
very  similar  to  those  presented  by  Hopensteadt  [16],  where  the  Cook 
model  was  used  with  an  added  constant  latent  period.  The  only  difference 
is  that  our  model  shows  the  successive  generations  of  cases  in  the 
beginning,  which  is  an  improvement. 
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CHAPTER  III 
MODEL  WITH  CONTROL 

3. I  Introduction 

In  his  presidential  address  to  the  Royal  Society  of 
Tropical  Medicine  and  Hygiene  (1965),  Dr.  G.  McDonald  [28],  while  talk¬ 
ing  of  mathematical  models  for  disease  spread,  said:  "When  it  comes  to 
development  of  policy  in  such  a  way  as  to  gain  the  initiative  and 
ultimately  gain  complete  mastery  over  an  infection  throughout  substantial 
areas  of  populations,  scientific  knowledge  has  as  yet  scarcely  been 
brought  into  play.  There  are  techniques  by  which  such  policy  can  be 
developed.  An  important  one  is  the  study  of  dynamics  of  infections,  with 
full  evaluation  of  their  relative  susceptibility,  in  terms  of  prevalence, 
to  alterations  in  different  aspect  of  environment.  In  this  way  their 
sensitivity  to  control  can  be  put  on  a  scientific  footing  and  much  of 
the  present  element  of  "Hit  and  Miss"  removed  from  the  general  study." 

The  remarks  quoted  above  underline  the  fact  that  any  advance 
in  the  art  of  mathematical  modeling  of  epidemics  is  in  itself  a 
contribution  to  the  control  of  disease  spread.  This  is  so  because  it 
helps  in  the  understanding  of  mechanism  of  disease  spread  and  hence 
helps  in  the  formulation  of  control  strategy.  In  this  respect  the  model 
discussed  in  the  previous  chapter  represents  a  significant  contribution, 
as  it  shows  the  sensitivity  of  the  disease  spread  to  initial  susceptibl es , 
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to  changes  in  3,  and  to  the  values  of  latent  and  infectious  periods. 

A  constant  monitoring  of  the  susceptible  population,  to  keep  it  below 
its  known  threshold  value,  and  control  of  3,  to  make  it  as  low  as 
possible,  go  a  long  way  to  prevent  the  disease  spread.  General 
notification,  closing  up  schools  and  other  public  places,  and  enforcement 
of  strict  quarantine  etc.,  are  the  methods  commonly  used  to  reduce  6 
and  to  prevent  disease  spread  when  such  a  danger  exists.  If  the 
quantitative  effect  of  these  measures  on  the  value  of  3  is  known,  their 
preventive  effect  can  be  easily  estimated  by  the  use  of  this  model. 

The  above  preventive  methods,  although  universally  used, 
have  only  a  limited  effect  and  are  difficult  to  enforce,  especially 
under  the  modern  conditions  of  swift  travel.  They  have  either  to  be 
supplemented  or  altogether  abandoned  in  favour  of  the  sophisticated 
methods  of  immunization  now  available.  These  methods,  which  create 
artificial  immunity  or  resistance  to  the  disease,  are  active  immunization, 
passive  immunization,  and  chemoprophylaxis.  This  chapter  is  primarily 
concerned  with  the  modification  of  the  epidemic  model  to  include  para¬ 
meters  representing  controls  by  the  first  two  of  these  methods.  An 
analysis  of  the  model  with  respect  to  the  sensitivity  of  its  response 
to  these  controls  is  presented. 

3.2  State  of  the  ccrt 

Although  a  considerable  amount  of  work  has  been  done  in  the 
field  of  mathematical  modeling  of  epidemics,  yet  the  incorporation  of 
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control  effects  in  these  models  is  only  in  the  beginning  stages.  In 
fact,  only  three  serious  efforts  where  control  is  incorporated  in  the 
epidemic  models  have  come  to  this  author's  notice;  those  of  Jacquette 
D7],  Taylor  [40]  and  ReVelle  [33].  Among  these,  the  last  is  the  one 
closest  to  the  situation  in  practice,  as  it  studies  the  quantitative 
effect  of  various  controls,  applied  at  different  rates,  on  the  final 
cost  of  the  disease  control.  The  other  two  studies  evaluate  the  economic 
benefits  of  stopping  a  disease  instantaneously  by  massive  application  of 
control.  The  assumption  of  unlimited  control  for  stopping  the  disease 
may  be  practical  from  the  point  of  view  of  disease  control  in  isolated 
cattle  herds,  but  is  hardly  realistic  from  the  public  health  administration 
point  of  view. 

In  practice,  the  control  of  disease  spread  in  large 
population,  always  strains  the  resources  of  public  health  authorities. 
Modern  practice  of  epidemic  control  consists  of  the  administration  of 
one  or  more  of  the  three  possible  controls:  (1)  active  immunization, 

(2)  passive  immunization  and  (3)  chemoprophylaxis  depending  upon  the 
disease  in  question.  Available  quantities  of  drugs  and  vaccines,  their 
rates  of  production  by  pharmaceutical  industry,  and  the  available 
number  of  public  health  personnel  to  administer  these  controls, are 
always  limited.  Hence  the  administration  of  these  controls  is  anything 
but  instantaneous.  Therefore,  in  the  interest  of  evolving  a  useful 
control  strategy,  it  is  imperative  that  the  sensitivity  of  the 
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epidemic  models  to  these  controls  is  studied  in  some  detail.  ReVelle's 
[33,  34,  23]  model,  though  quite  realistic,  is  not  general  enough  (it 
was  developed  specifically  for  tuberculosis)  and  has  not  been  used  for 
the  study  of  sensitivity  of  the  disease  to  the  controls  (B.C.G. 
vaccination  and  chemoprophylaxis  in  the  case  of  tuberculosis).  To 
remedy  this  situation,  the  model  presented  in  Chapter  II  is  modified, 
in  the  next  section,  to  include  active  and  passive  immunization  controls, 
and  its  sensitivity  to  each  of  the  two  controls  is  analysed.  The  model 
thus  obtained  is  similar  to  the  one  presented  previously  by  Gupta  and 
Rink  [13]. 

3.3.  Active  Immunization 

Active  immunization  is  the  process  of  creating  artificial 
immunity  in  the  system  of  a  susceptible  by  injecting  a  vaccine  of  either 
dead  or  live,  but  attenuated,  disease  micro-organi sm.  The  method 
induces  the  body's  defense  mechanism  to  produce  antibodies  specific  to 
the  micro-organism  of  the  disease.  The  experience  of  producing  the 
antibodies,  thus  gained  by  the  susceptible,  confers  on  him  immunity  to 
the  disease  in  question.  This  immunity  is  almost  comparable  to  the  one 
gained  on  recovery  from  the  disease.  Its  duration  depends  on  the 
disease  in  question  and  the  method  used  for  the  preparation  and  admi¬ 
nistration  of  the  vaccine. 

The  effect  of  administering  active  vaccine  to  susceptibles 
is  that  of  transferring  them  from  the  susceptible  population  to  the 
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immune  population,  directly,  without  their  going  through  the  cycle  of 
disease.  An  ideal  situation  would  be  to  immunize  every  susceptible, 
thus  eliminating  the  danger  of  disease  spread  altogether.  This  ideal 
is  practically  impossible  to  achieve  and  is  economically  undesirable, 
as  almost  the  same  results  can  be  obtained  by  a  constant  surveillance 
of  the  population  and  administration  of  vaccine  at  such  a  rate  that  the 
population  of  susceptibles  is  always  kept  below  the  threshold  value  for 
the  disease.  Considering,  however,  the  wide  variety  of  diseases  that 
may  strike  a  population  and  the  general  reluctance  or  negligence  on  the 
part  of  individuals  to  get  immunization  when  the  danger  is  not  imminent, 
it  is  almost  impossible  to  have  sustained  active  immunization.  There¬ 
fore,  there  is  always  a  possibility  of  an  epidemic  outbreak  catching  a 
community  inadequately  protected,  and  control  has  to  be  applied  after 
the  first  signs  of  epidemic  appear.  It  is  this  situation  for  which  the 
model  is  modified  to  account. 

One  way  of  representing  the  active  control  is  to  assume 
that  vaccine  is  given  only  to  known  susceptibles  (which  can  be  done  by 
giving  a  test  before  giving  the  actual  vaccine,  as  in  the  case  of  R.C.G.) 
and  each  dose  is  effective  in  providing  immunity,  immediately  after 
administration.  Thus  if  U-|(t)  is  the  rate  at  which  susceptibles  are 
being  vaccinated  at  time  t  the  susceptible  population  can  be  said  to 
be  reducing  at  that  rate,  due  to  active  control.  This  fact  can  be  in¬ 
corporated  in  the  model  by  changing  equation  (2.4)  from 
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X1  ‘  6o  X1  X2  (3.1 ) 

t0  X1  =  -  B0  X1  X2  -  Ul(t)-  (3.2) 

The  assumption  that  an  active  vaccine  is  effective 
immediately  after  its  administration  is  not  wholly  realistic,  as  all 
active  vaccines  need  some  latent  period  of  their  own  for  the  immunity 
response.  This  may,  however,  be  a  reasonable  assumption  in  the  case  of 
those  diseases  which  have  a  long  latent  period,  since  in  such  cases 
immunity  can  develop  before  the  disease  does,  even  with  a  concurrent 
exposure. 

Correspondingly ,  in  the  case  of  the  normalized  model 
equation  (2.17)  is  changed  from 


*1  =  -3  x1  x2 

(3.3) 

to 

x i  —  “6  x i  x 2  -  u i  ( t ) , 

(3.4) 

where  u-^  =  U-j  is  the  number  of  vaccinations  given  per  day,  expressed 
as  a  fraction  of  the  total  population,  and  is  the  per  unit  equivalent 
of  U-| . 

Since  it  is  not  always  convenient  and  even  in  some  cases  not 
possible  to  give  a  test  for  susceptibility,  a  more  general  situation  will 
be  where  vaccine  is  given  at  random.  In  such  a  case  we  can  assume  that 
of  the  U-|(t)  vaccines  given  per  day,  at  time  t,  only  the  fraction  proportio¬ 
nal  to  the  susceptible  population  is  useful  and  the  others  are  wasted. 

So  the  effective  immunization  rate  would,  in  that  case,  be  U-|(t).  X-j 

N“ 
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for  un-normal  ized  case  and  u-|(t)  x-j  for  normalized  case. 

It  has  been  assumed  so  far  that  every  vaccination  given  to 
a  susceptible  is  a  success.  This  is  not  true  in  actual  practice,  as  we 
know  now  that  there  is  significant  failure  rate  which  differs  from 
vaccine  to  vaccine.  This  situation,  however,  can  be  easily  considered 
by  assuming  U-|(t)  as  the  number  of  effective  vaccinations  rather  than 
actual  vaccinations  per  day.  Another  possibility  to  be  considered  is 
the  delay  in  the  effectiveness  of  vaccine.  If  an  active  vaccine  may, 
on  the  average,  be  considered  effective  tc  days  after  it  is  administered, 
the  control  term  becomes  U]  (t~rc)  or  U]  (t-r  ),  as  the  case  may  be. 

The  equations(3.2)  and  (3.4)  respectively ,  become 

•  X 

xi  =  -  Bi  xi  x2  -  ut  •  r  0.5) 

x-j  =  -3  x i  x2  -  u1  ( t — t c )  •  X-, .  (3.6) 

These  equations,  when  incorporated  in  the  respective  versions  of  the 
model,  give  the  modified  control  model  for  active  control.  As  an 
example,  the  normalized  model  with  active  control  will  be  as  follows: 


*1 

3  X-j  x^  -  X-j  u-j  (t-x^) 

(3.7) 

x2 

=  aft)  -  r(t) 

(3.8) 

X3 

=  r(t) 

(3.9) 

x4 

=  3  K(t)  •  X-j  x2 

(3.10) 

x5 

=  u1 (t) 

(3.11) 
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where 


a(t)  = 


1 


rr'  ~  i  *4<t'T) 

/  2„  Jq 


2u 


1 


exp 


2a- 


(y-j-x) 


r(t)  = 


1 


rr-'  A  I  ^(t'x) 

■rs  a2  j 


2u. 


exp 


2  a 


1  ( 

— 2  i^-xj 


dx 

(3.12) 


dx , 
(3.13) 


and  other  forms  of  a(t)  and  r(t)  are  the  same  as  those  presented  in 
Chapter  II,  equations  (2.23)  to  (2.26).  The  new  variable  x,-(t)  repre¬ 
sents  the  fraction  of  population  actively  immunized  up  to  time  t; 
x^(t^)  will,  therefore,  represent  the  total  vaccine  used  at  final 
time  t^. 


Results  of  the  solution 


The  model  is  now  analysed  for  different  values  of  active 
control  u -j  •  Control  delay  xc  can  be  assumed  to  be  equal  to  zero  with¬ 
out  loss  of  generality,  as  the  control  is  first  applied  on  a  certain  day 
and  any  delay  xc  simply  shifts  this  time  of  application. 

Figure  3.1  shows  the  effect  of  variable  control  u-|  applied 
over  a  fixed  interval  of  time.  Administration  of  control,  at  rates 
shown  on  each  curve,  is  assumed  to  start  on  the  15th  day  (about  one  mean 
incubation  period  after  the  start  of  disease)  and  end  on  the  25th  day. 

We  observe  from  the  figure  that  the  higher  the  rate  of  active  control  , 
the  lower  the  total  number  of  cases. 


Figure  3.2  shows  the  effect  of  different  control  rates  with 
the  total  amount  of  control  administered  remaining  the  same  in  each  case 
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FIG.  3*2  EFFECT  OF  CONSTANT  ACTIVE  CONTROL  APPLIED 
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and,  again,  the  control  starting  at  the  15th  day,  but  this  time  ending 
on  a  different  day.  The  plot  of  the  disease  response,  in  this  case, 
shows  that  best  results  are  obtained  by  using  the  control  at  the  high¬ 
est  rate  although  the  total  amount  of  control  remains  the  same. 

Figure  3.3  shows  the  effect  of  delay  in  the  application  of 
control.  In  each  case  the  same  amount  of  control  (0.3)  is  used  and  at 
the  same  rate  (u-j  =  0.03),  but  the  control  starts  at  different  times. 
The  delay  in  the  start  of  the  administration  of  control  represents  the 
time  lag  in  putting  the  control  effort  in  action.  This  time  lag  may 
be  due  to  the  immediate  non-availability  of  the  vaccine  in  sufficient 
quantities,  or  simply  due  to  a  failure  to  realize  the  epidemic  nature 
of  the  disease.  We  find  from  figure  3.3  that  the  disease  is  best 
controlled  by  the  earliest  action,  and  the  larger  the  delay,  the  less 
effective  is  the  control. 

Thus  we  conclude  from  these  three  plots  that  the  best 
results  are  obtained  by  applying  the  available  active  control  at  the 
highest  possible  rate  and  at  the  earliest  available  opportunity.  These 
results  seem  to  be  what  common  sense  would  suggest.  However,  since  now 
we  are  able  to  evolve  a  quantitative  measure  of  the  effect  of  active 
control,  the  feasibility  of  an  optimum  control  strategy  is  very  clear, 
and  the  next  chapter  evolves  this  strategy  by  using  Pontryagin's 
Minimum  Principle. 


' 


DISEASE  REPORT  RATE  (*10‘ 


-49- 


FIG.  3-3  EFFECT  OF  DELAY  IN  APPLICATION  OF 

ACTIVE  CONTROL 


-  50  - 


3.4  Passive  Immunization 

As  opposed  to  active  immunity,  which  is  the  process  of 
creating  antibodies  in  the  system  of  a  susceptible,  passive  immunity 
is  provided  by  injecting  "ready  made"  antibodies  into  the  system  of  an 
individual.  This  form  of  immunity  is  immediately  effective  and  is  even 
curative  in  nature.  The  immunity  provided  in  this  manner,  however,  is 
not  as  permanent  as  one  conferred  by  active  immunization,  and  is  liable 
to  be  lost  quite  soon.  Other  hazards  of  passive  immunization  are  the 
possibility  of  introduction  of  other  infections  from  the  source  of 
antibodies  (human  or  animal  blood  plasma)  and  incompatibility  of  the 
system  to  foreign  plasma.  Notwithstanding  the  above  difficulties, 
passive  immunization  is  now  an  accepted  method  of  control  for  those 
diseases  for  which  antiserum,  needed  for  passive  immunization,  can  be 
easily  produced.  This  antiserum  is  generally  administered  to  those 
suspected  of  already  carrying  the  disease  micro-organism  in  their  system, 
thus  making  them  ineffective  in  spreading  the  disease  further. 

The  accurate  quantitative  effect  of  the  passive  immunization 
on  the  latent  and  infectious  periods  of  the  disease,  and  the  duration  of 
immunity  it  may  confer  on  those  immunized,  is  not  yet  fully  known.  The 
methods  of  representation  of  passive  control  discussed  in  the  next  few 
paragraphs  may,  therefore,  not  be  the  best  possible,  but  may  serve  as 
a  guide.  More  over  these  methods  of  representation  may,  hopefully,  also 
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be  used  to  represent  the  action  of  chemoprophyl axis  a  method  of 
control  which  is  not  discussed  in  this  thesis  for  lack  of  knowledge 
about  its  action  on  various  sub-populations  of  the  model.  Two 
different  methods  of  representation  of  passive  immunization  have  been 
used  to  adapt  the  model  to  passive  control.  Sensitivity  of  the  model 
to  changes  in  the  passive  control  has  been  analysed  in  both  cases. 

Case  t 

In  this  case  the  passive  immunization  is  assumed  to  be 
administered  at  random,  thus  removing  both  susceDtibles  and  infectives 
from  circulation  by  conferring  on  them  a  temporary  immunity.  Represent¬ 
ing  the  passive  control  by  U2  inoculations  per  day,  in  the  same  way  as 
in  the  case  of  the  active  control  model  but  with  a  difference  that  now 
immunity  is  conferred  both  to  susceptibles  and  infectives,  the  modified 
model  for  the  normalized  case  becomes 


X1 

=  -  6  x-j  -  u^  x-j 

(3.14) 

x2 

=  a ( t )  -  r(t)  -  u2  x2 

(3.15) 

x3 

=  r(t) 

(3.16) 

X4 

=  K(t)  -8x^2 

(3.17) 

X6 

=  u2 

(3.18) 

where  u2  =  and  are  the  same  as  those  given  by  equations 
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(3.12)  and  (3.13)  respectively,  and  is  the  new  variable  representing 
the  total  amount  of  passive  control.  It  may  be  noted,  here,  that  no 
time  delay  is  considered  in  the  control  variable  in  this  case,  since 
it  is  known  that  passive  control  acts  almost  instantaneously.  In  fact 
this  is  the  major  advantage  of  using  this  control. 

Results  of  the  solution 

Figure  3.4  represents  the  solution  of  the  model  for  differ¬ 
ent  rates  of  passive  control.  We  find  that  the  higher  the  control  rate 
used,  the  smaller  the  number  of  cases.  It  may  also  be  noted  that,  for 
the  same  reduction  in  cases,  the  amount  of  u^  used  is  much  smaller  than 
the  amount  of  u^.  This  is  an  expected  result,  as  we  know  that  passive 
control  acts  on  the  infective  cases  also, and  as  such  is  a  more  effective 
way  of  suppressing  disease,  at  least  temporarily. 

Figure  3.5  shows  the  effect  of  delay  in  application  of 
passive  control.  A  constant  (.01)  rate  of  passive  control  is  used  in 
all  the  cases  but  it  is  started  at  different  delays  after  the  start  of 
disease.  We  note  that  the  effect  of  control  is  maximum  when  it  is 
applied  at  the  earliest  time. 

Case  2 

In  case  1  it  was  assumed  that  the  passive  immunization  is 
given  to  the  population  at  random.  However,  in  practice  it  is  prefer¬ 
able  to  give  passive  immunization  to  the  known  contacts  of  reported 
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FIG.  3-4  EFFECT  OF  VARIABLE  PASSIVE  CONTROL 

CASE  1 
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EFFECT  OF  DELAY  IN  PASSIVE  CONTROL  APPLICATION 

CASE  1 


FIG.  3-5 
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cases.  Thus,  to  represent  a  situation  of  this  kind  we  assume  that 
when  a  case  is  reported  its  contacts  are  traced  and  those  who  might 
have  picked  up  the  disease  may  be  given  passive  immunization  to 
suppress  the  disease.  We  can  assume  that  on  the  average  U2  infectives 
are  given  passive  immunization,  per  reported  case.  Thus  the  rate  of 
reduction  in  the  infective  population  due  to  passive  vaccination  may 
be  assumed  to  be  R(t)  per  day.  This  rate  becomes  U2  r(t)  for  a 
normalized  model.  It  may  be  noted  here  that  in  this  case  is  a  non- 
dimensional  quantity,  as  it  represents  the  number  cf  infectives  traced 
per  reported  case. 

The  modified  normalized  model  now  is: 


X1 

:  -  M|  x2 

(3.19) 

x2 

=  aft)  -  r(t)  -  U2  r(t) 

(3.20) 

X3 

=  r(t) 

(3.21) 

X4 

=  K(t)  6  X-|  x 2 

(3.22) 

X6 

=  U2  r(t) , 

(3.23) 

where  a(t)  and  r(t)  are  as  defined  in  case  1  and  xc  corresponds  to  xc. 

o  b 

Results  of  the  solution 

Figure  3.6  shows  the  result  of  the  solution  of  the  model  for 
various  values  of  U2-  We  find  that  as  U^,  the  number  of  infectives  re¬ 
moved  per  reported  case,  is  increased  the  disease  is  better  controlled. 
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FIG.  3-6  EFFECT  OF  VARIABLE  PASSIVE  CONTROL 

CASE  2 
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Figure  3.7  shows  the  effect  of  delay  in  starting  the  admi¬ 
nistration  of  passive  control.  As  would  be  expected,  the  control  is 
most  effective  when  started  at  the  earliest  time.  Thus  we  conclude, 
again,  that  for  best  results  maximum  control  effort  should  be  applied 
at  the  earliest  possible  time. 

3.  5  Conotus'ions 

In  this  Chapter  the  sensitivity  of  the  model  to  active  and 
passive  control s , each  taken  separately,  has  been  analysed.  Control  with 
chemoprophylaxis  was  not  considered.  It  is,  however,  clear  that  the 
third  control  can  also  be  incorporated  in  the  same  manner  as  the  other 
two.  While  considering  the  model  with  control,  no  limitations  on  the 
availability  of  the  vaccines  and  antisera,  and  on  their  costs,  have  thus 
far  been  considered.  In  practice,  however,  the  more  important  question 
is:  When  more  than  one  control  method  is  available,  how  much  of  each 
should  be  used  for  best  results  ?  This  problem  is  very  difficult  to 
solve  with  the  methods  of  analysis  used  in  this  chapter,  since,  if 
various  combinations  of  controls  were  to  be  tried  in  a  hope  of  finding 
the  best  combination,  we  are  likely  to  end  up  doing  a  large  number  of 
solutions  of  the  model,  yet  having  no  certainty  of  finding  the  best  com¬ 
bination  of  controls.  This  problem  will,  however,  be  solved  in  the  next 
chapter  using  the  method  of  dynamic  optimization. 


In  conclusion,  as  a  result  of  the  discussion  so  far,  we  can 
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FIG.  3-7  EFFECT  OF  DELAY  IN  PASSIVE  CONTROL  APPLICATION 

CASE  2 
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write  the  model  incorporating  both  the  controls.  Since  we  have  used 
two  different  methods  of  representing  passive  control,  we  have  two 
different  forms  of  the  model.  Only  the  normalized  model  is  given  here 
for  both  the  cases. 


Case  l 


Case  2 


x-,  =  ■  0  x1  x2  -  ui(t-'rc)  •  X1  ■  u2  X1  (3.24) 

x2  =  a(t)  -  r(t)  -  u2  x2  (3.25) 

x3  =  r(t)  (3.26) 

x4  =  K(t)  -0x^2  (3.27) 

x5  =  u1 (t)  (3.28) 

x6  =  u2(t).  (3.29) 

xT  =  -  6  x1  x2  -  u]  ( t~T c )  •  X-|  (3.30) 

x 2  =  a(t)  -  r (t)  -  U2  r (t)  (3.31 ) 

x3  =  r (t)  (3.32) 

x4  =  K(t )  •  6  x-j  x2  (3.33) 

x5  =  u-j(t)  (3.34) 

x 6  =  U2  r (t)  .  (3.35) 
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In  both  the  cases  a(t)  and  r(t)  are  defined  as  in  equations 
(2.21)  to  (2.26)  in  Chapter  II,  depending  on  the  choice  of  represent¬ 
ation  to  be  used.  Representative  expressions  for  a(t)  and  r(t)  when 
mean  value  of  incubation  period  is  used  are: 


a(t)  =  -L  .  1 

/  2tt  a 


2p, 


1  -0 


x,  (t-T)  •  exp 


1  (yi-f)2 


2  a 


2  VM1 


1 


dt  (3.36) 


r(t)  = 


1  1 


/  2t r 


a. 


2 

x4  (t-T )  •  exp 

r  i  /  x2l 

-  2  (v2-t) 

J 

0 

2“2  -1 

dr.  (3.37) 
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CHAPTER  IV 

OPTIMUM  CONTROL  MODEL 

4. I  Introduction 

Constant  vigil  kept  by  the  public  health  organizations 
around  the  world  has  helped  to  allay,  considerably,  the  danger  of 
full-scale  epidemic  outbreaks.  This  preventive  effort,  however,  entails 
large  sums  of  money,  and  efforts  are  always  underway  to  improve  the 
control  procedures  so  as  to  reduce  this  expenditure  and  get  maximum 
benefits  with  minimum  outlay.  The  mathematical  theory  of  epidemics, 
by  providing  a  better  insight  into  the  mechanisms  of  disease  spread, 
has  indirectly  contributed  to  the  fulfilment  of  the  above  stated  goal, 
since  a  better  understanding  of  the  mechanisms  of  disease  spread* 
always  ,resul ts  in  improved  methods  for  its  prevention.  There  is,  how¬ 
ever,  a  consensus  in  the  literature  that  the  mathematical  models  have 
not  yet  been  used  to  their  fullest  potential  for  the  control  of  epi¬ 
demics.  In  particular,  the  use  of  these  models  in  determining  a  most 
economical  control  strategy  between  various  competing  controls,  by  the 
use  of  dynamic  optimization  techniques,  has  been  frequently  predicted. 
This  chapter  presents  a  procedure  for  the  application  of  dynamic 
optimization  theory  to  the  control  model  developed  in  Chapter  III. 

Three  references,  viz.  ReVelle  [33],  Taylor  [40],  and 
Jacquette  [17],  were  discussed  in  the  last  chapter  for  their  contri- 
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bution  to  the  field  of  epidemic  control  models.  It  is  significant  to 
note  that  the  same  three  research  efforts  are  also  the  ones  that  con¬ 
sider  the  question  of  economical  control  strategy,  and  no  other  work 
in  this  field  has  yet  come  to  the  notice  of  this  author.  ReVelle 
[33]  considered  three  competing  controls  (prophylaxis,  cures  and 
B.C.G.  vaccination)  for  his  tuberculosis  model.  Four  alternative 
control  strategies,  each  reducing  the  epidemic  to  an  assumed  level  in 
a  fixed  time,  were  considered,  and  the  one  involving  minimum  control 
cost  was  identified.  Taylor  [40]  considered  a  herd  of  dairy  cattle 
which  can  be  immunized,  at  any  time,  to  confer  complete  immunity,  if 
so  desired.  A  vaccination  schedule  that  minimizes  the  long-run  time- 
average  sum  of  the  costs  of  immunization  and  expected  disease  losses 
was  determined.  Jacquette  H7],  in  turn,  assumed  that  an  epidemic  can 
be  stopped  instantly  by  a  massive  dose  of  immunization  and  proceeded 
to  find  out  the  most  economical  time  of  this  stopping. 

Although  the  three  studies,  summarized  above,  are  signifi¬ 
cant  contributions  to  the  field  of  economical  control  of  epidemics, 
yet  each  stopped  short  of  the  stated  ideal  of  determining  control 
strategy  by  dynamic  optimization  techniques.  The  following  quotation 
from  the  latest  of  these  studies  [17,  pp.  11-12]  illustrates  the  point 
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"We  will  see  that  the  standard  calculus  of  variation  or 
control  theory  techniques  for  dynamic  control  of 
coefficients  (or  state  variables)  fail  since  the 
differential  equations  of  the  model  are  second  degree. 
Ideally  a  solution  consisting  of  an  optimal  control 
trajectory  is  desired,  but  this  must  be  left  as  a 
topic  for  further  work.  Steady  state  solutions  are 
available  and  controllable  coefficients  can  be  set  to 
minimize  the  averaqe  cost  per  unit  time.  In  cases 
where  the  model  variable  is  a  known  function  of  time 
we  will  look  at  a  direct  form  of  control  and  discover 
an  analogy  with  inventory  theory." 

The  new  model,  discussed  in  the  previous  chapters,  removes 
some  of  these  difficulties  and  is  used  for  determining  an  optimum 
control  strategy  by  the  direct  application  of  Pontryagin's  Minimum 
Principle.  The  model  has  been  converted  to  a  set  of  non-linear 
differential  difference  equations  by  the  substitution  of  the  variables 
a ( t)  and  r(t)  in  their  discrete  form  (equations  2.24  and  2.25).  The 
resulting  state  equations  are  functions  of  delayed  states  and  delayed 
control,  both  with  finite  number  of  delays.  Solutions  resulting  from 
the  optimum  control  study  of  this  model  are  presented. 


4.2  Formulation  of  the  control  problem 


The  final  form  of  the  normalized  control  model,  incorpo¬ 
rating  both  active  and  passive  immunization  controls,  was  presented  in 
the  concluding  section  of  Chapter  III.  This  model  can  now  be  modified, 
by  reducing  its  dimensionality,  to  make  it  more  suitable  for  optimiz¬ 
ation  study.  It  may  be  pointed  out  here  that  the  results  obtained 
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from  the  normalized  model,  which  will  be  the  only  form  considered  here, 
can  easily  be  extended  to  non-normal ized  model,  if  so  desired. 

A  quick  review  of  equations  (3.24)  to  (3.37),  representing 
the  model  for  two  different  formulations  of  passive  control,  reveals 
that  the  fourth  equation  of  the  model,  in  both  the  cases,  can  be  easily 
eliminated  by  substituting  it  in  the  expressions  for  variables  a(t)  and 
r(t).  Moreover,  by  introducing  a  new  equation  representing  the  cost  of 
the  control  used,  the  last  two  equations  of  the  model  can  also  be 
dropped  without,  in  any  way,  reducing  the  effectiveness  of  the  model 
for  optimization.  Since  it  is  intended  to  obtain  the  model  in  the 
differential -difference  form  rather  than  the  integro-differential  form, 
discrete  versions  of  a(t)  and  r(t)  (represented  by  equations  (2.24)  and 
(2.25)  rather  than  those  represented  by  equations  (3.36)  and  (3.37)), 
will  be  used  in  the  modified  model.  In  addition,  we  can  replace  the 
' x 1  variables  by  y  to  avoid  confusion  in  the  subscripts  used  in  the 
previous  chapters. 

Thus  keeping  the  above  modifications  in  mind  and  substi¬ 
tuting  the  resulting  expressions  for  a(t)  and  r(t)  in  the  main  equations 
of  the  model,  the  model,  for  two  cases  of  the  last  chanter,  now  becomes: 

Case  l 
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y4  =  A  u-|(t)  +  B  u-j2(t)  +  C  u2(t)  +  D  u22(t) 


(4.4) 
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y4  =  A  u]  (t)  +  B  u^  (t)  +  {C  +  D  U2(t)}  U2(t)  x 

2y 

Ke  1  \  r  r  1  /  ,2, 

2^  exp  { - 2  (y2-nr)  } 


/  2tt 


u2  n=o  L  c02 

y1  (t-m)  •  y2(t-nT)  . 


(4.8) 


Where  y-j  ,  y^  and  y^  represent  the  per  unit  sub-populations  of 
susceptibl es ,  infectives  and  recovered  or  removed  cases;  y4  is  the 
cost  of  immunization  expressed  on  per  unit  basis;  xc  is  the  delay,  in 
days,  in  the  action  of  active  control,  and  t  is  the  unit  of  time,  one 
day  in  this  case.  Again,  A  and  C  are  the  linear  and  B  and  D  the 
quadrature  costs  for  the  active  and  passive  rates  of  immunization, 
respectively.  A  more  explicit  discussion  of  constants  A,  B,  C  and  D 
will  be  given  in  a  later  section,  where  the  numerical  values  for  these 
constants  are  chosen. 


Using  the  optimum  control  terminology,  we  can  now  call 
y i»  y2>  y3  and  y4  the  states,  and  equations  (4.1)  to  (4.4)  (also 
equations  4.5  to  4.8)  the  state  equations  of  the  model.  If  the  initial 
values  of  states  and  their  history  up  to  the  initial  time  t  is  known, 
the  state  equations  can  be  solved  for  a  given  control  function,  using 
the  procedure  described  in  Chapters  II  &  III.  At  the  final  time  t*, 
we  then  obtain  the  final  values  of  the  states,  i.e.  susceptibl es , 
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infectives,  recovered  or  removed  cases,  and  the  cost  of  immunization 
used  during  the  course  of  the  epidemic. 

Cost  function 

The  final  cost  of  an  epidemic  can  be  assessed  by  adding 
the  cost  of  cases  afflicted  to  the  cost  of  control  used  during  the 
epidemic  period.  It  is  not  suggested  here  that  cost  to  an  individual 
struck  with  disease  can  be  counted  in  dollars  and  cents,  especially 
in  a  fatal  case,  but  what  we  mean  by  cost  of  a  case  is  the  weighted 
average  cost  with  respect  to  the  control  cost  (which  can  be  determined 
more  accurately).  In  physical  terms,  this  weighted  average  cost  may 
represent  the  cost  of  average  man-hours  lost  plus  the  cost  of  hospital¬ 
ization,  for  each  reported  case.  In  addition,  it  may  be  argued  that 
if  the  disease  is  well  controlled  the  fatality  rate  will  be  negligibly 
small  and  hence  the  cost  per  case,  calculated  on  the  basis  defined 
above,  is  more  near  to  the  true  cost.  When  the  model  is  applied  to  an 
epidemic  in  a  herd  of  cattle,  however,  the  cost  per  case  is  more 
straight  forward.  Let  us,  therefore,  assume  that  the  cost  of  one 
reported  case  is  C-j .  The  cost  function,  on  the  per  unit  population 
basis  can  then  be  written  as 


J  =  C1  y3(tf)  +  y4(tf) 


(4.9) 
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where  .3  is  the  total  cost  of  the  controlled  epidemic;  y^(t^)  is  the 
final  value  of  the  reported  cases,  and  y4(tf)  is  the  total  cost  of 
the  control  applied  during  the  course  of  the  epidemic.  Equation  (4.9) 
can  be  modified  to  also  include  the  cost  for  the  residual  infectives, 
if  any,  at  the  final  time,  thus  giving, 

J  =  C1  y3(tf)  +  C1  y2(tf)  +  y4(tf).  (4.10) 

A  control  policy  which  minimizes  the  cost  J  at  the  end  of 
an  epidemic  is  the  optimum  policy.  This  control  policy  can  now  be 
formulated  by  the  application  of  Pontryagin's  Minimum  Principle. 

We  thus  have  a  vector  state  equation  which  consists  of  a 
set  of  non-linear  differential -difference  equations  with  a  finite 
number  of  delays  in  both  state  and  control.  The  cost  function  is  a 
function  of  the  final  values  of  the  state  only,  thus  giving  a  Meyer 
problem.  The  optimization  problem,  although  difficult,  yet  is  not 
impossible  to  solve,  especially  in  view  of  the  latest  research  reported 
in  the  field.  A  closed-form  solution  is,  however,  impossible  at  this 
stage,  and  hence  only  a  numerical  solution  is  attempted. 

4 .  3  State  of  the  art 

A  considerable  amount  of  work  has  already  been  done  in  the 
field  of  optimization  of  time  delay  control  systems.  Yet  the  theory 
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relating  to  this  subject  is  being  continuously  updated  as  the  solutions 
to  an  increasing  number  of  problems  of  this  nature  are  attempted.  This 
fact  is  evidenced  by  some  recent  references  in  the  field:  [36],  [27], 
[9]. 

Direct  application  of  Pontryagin's  Maximum  (or  Minimum) 
Principle  to  systems  with  time  delay  is  due  to  G.L.  Kharatishvil i . 

He  first  applied  the  Maximum  Principle  to  a  system  having  a  delay  in 
state  but  not  control;  see  Pontryagin  et  al  [32,  p.  213].  Subsequent¬ 
ly,  the  method  was  extended  by  him  to  include  the  case  of  a  delay  both 
in  state  and  in  control;  see  Kharatishvil i  [22],  More  recently  (1970), 
Budelis  and  Bryson  [9]  developed  necessary  conditions  for  an  extremal 
path  in  the  case  where  the  system  equations  and  the  performance  index 
contain  a  time  delay  in  both  the  state  and  the  control  variables.  An 
analytical  solution  for  a  linear  system  with  a  quadratic  performance 
index,  where  control  variable  appears  in  the  system  equations  at  the 
present  time  and  at  a  previous  time,  was  also  presented  by  Budelis  and 
Bryson  [9].  McAulay  [27]  derived  the  optimum  control  criterion  for  a 
system  consisting  of  non-linear  differential -difference  equations  with 
a  finite  number  of  delays  in  the  state.  He,  furthermore,  extended  the 
application  of  steepest  descent  (gradient)  technique  of  Kelly  [20]  and 
Bryson  [7]  to  this  case. 
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Concerning  the  gradient  methods,  Gottlieb  [12]  has  brought 
out  the  relationship  of  these  methods  with  the  calculus  of  variations 
and  Min-H  (Pontryagin ' s  Minimum  Principle)  methods.  He  has  also 
proposed  two  rapid  convergence  gradient  methods  for  solving  Min-H 
problems.  Lasdon  et  al  [25]  have  extended  the  application  of  con¬ 
jugate  gradient  method  to  the  Min-H  problem,  thus  proving  the  applica¬ 
bility  of  general  parallel  tangent  methods  to  this  problem.  Shah  et 
al  [38]  have  applied  a  parallel  tangent  (abbreviated  "PARTAN")  method 
to  the  minimization  of  a  function  of  several  variables. 

A  reasonably  comprehensive  solution  of  the  epidemic  control 
problem,  formulated  in  the  last  section,  can  be  obtained  using  the 
information  contained  in  the  research  reviewed  above.  Although  some  of 
the  theories  developed  in  the  research  detailed  above  have  been  success¬ 
fully  applied,  yet  many  of  them  have  not,  so  far,  been  rigorously 
tested  with  a  practical  application.  As  such,  the  present  problem 
provides  an  excellent  practical  application  of  the  techniques  discussed 
above. 

4.4  Optimization  procedure 

The  simulated  epidemic  control  model  was  solved  in  Chapter 
III  with  known  initial  conditions  and  a  termination  criterion  depending 
on  the  final  values  of  some  states.  Thus,  the  optimum  control  problem. 


P  l'UT  ) 


formulated  in  section  4.2,  is  essentially  one  of  fixed  initial  time 
and  free  terminal  time.  In  view  of  the  mathematical  complication 
introduced  by  the  free  final  time,  the  problem  is,  here,  simplified 
by  converting  it  to  one  of  fixed  initial  and  fixed  terminal  time. 

This  is  achieved  by  introducing  a  penalty  function  in  the  state 
equation.  The  procedure  applied  is  explained  in  the  following  paragraph. 

When  solving  the  optimum  control  problem  with  fixed 
terminal  time,  it  is  likely  that  absolute  minimum  of  the  cost  function 
would  be  obtained  with  the  final  values  of  some  states  becoming 
negative.  This  is  a  highly  undesirable  situation,  since  negative  sub¬ 
populations  have  no  physical  meaning.  In  other  words,  we  can  say  that 
there  is  an  implicit  constraint  on  the  state  variables,  and  this 
constraint  is  that  the  states  should,  always,  have  values  greater  than 
zero.  This  difficulty  can  be  overcome  by  using  a  penalty  constraint 
of  the  type  suggested  by  Kelly  [20,  p.  215].  Therefore  the  penalty 
function  needed  to  be  added  to  the  state  equation  for  y^  is 

1  Pj  yj2,  (4-1!) 

where  p.  is  a  positive  constant,  suitably  selected  for  the  correspond- 

\J 

ing  state,  and  6  is  a  Heaviside  unit  step  function  of  argument  y..  In 

vJ 

practice,  however,  we  find  that  for  the  present  problem,  if  a  heavy 
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penalty  constant  is  used  for  y2  (the  state  representing  the  infective 
population),  penalties  on  the  other  two  states  are  unnecessary. 

Concentrating  our  attention  on  case  1  of  the  model,  and 
incorporating  the  inequality  constraint  on  state  y^,  the  state 
equations  can  be  rewritten  as: 


Yi  = 
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y4  =  P2  6(y2)  •  y22  +  A  u1  (t)  +  B  u^(t)  +  C  u2(t)  +  D  u22(t) .  (4.1 5) 

The  above  state  equations  can  now  be  represented  in  the 
vector  form  as: 


y(t-nx),  u,  u(t-xc)j 


y  =  f [y.  y(t-x) ,  y(t-2x) ,  . . . 


(4.16) 
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where  y  is  a  state  vector  with  four  components  y-j  ,  y^,  y^  and  y^; 
f  is  the  corresponding  vector  of  functions  f^  ,  f ^ and  f^,  each 
representing  the  right  hand  side  of  one  of  the  above  four  equations; 
u  is  the  control  vector,  and  n  is  an  integer  equal  to  2^. 

From  physical  considerations  we  know  that  u(t)  is  piece- 
wise  continuous  with  finite  discontinuities  and  that  y(t)  is  con¬ 
tinuous  with  piecewise  continuous  first  derivatives. 


From  equation  (4.10)  we  see  that  the  cost  function  J  is  a 
function  of  the  final  values  of  states  at  preselected  final  time  t^. 
Mathematical ly 


J 


4>(y |  f). 


(4.17) 


thus  giving  a  Meyer  problem  which  satisfies  all  the  conditions  for 
application  of  PMP  to  systems  with  delay,  as  laid  down  by  Pontryagin 
et  al  [32],  Kharatishvil i  [22],  and  Budelis  et  al .  [9].  The  results 
of  Budelis  et  al  [9]  for  the  gradient  of  Hamiltonian  can  thus  be 
used  directly,  with  extension  to  a  multi-delay  system.  This  extension 
can  be  carried  out  on  the  basis  of  results  obtained  by  McAulay  [27], 
since  the  delays  encountered  in  equation  (4.16)  are  similar  to  those 
assumed  by  him. 


Therefore,  assuming  a  co-state  vector 
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X  [X-j,X2>  X3>X4],  (4.18) 

the  Hamiltonian  H,  for  the  above  system,  becomes 

H  =  X1f1  +  X2f2  +  X3f3  +  x4f4.  (4.19) 

The  costate  equations  of  the  system,  satisfying  the  con¬ 
ditions  of  optimality,  and  based  on  the  results  of  references  [9] 
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and 


(4.21) 


Moreover,  the  necessary  conditions  for  an  extremal 
(dJ=o,  for  arbitrary  6u)  are  that  H(t)  be  continuous  and  that 
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(4.22  a) 


tf“Tc  <  t  -  tf 


(4.22  b) 


The  problem  is  now  solved  by  the  use  of  parallel  tangent 
gradient  method.  The  numerical  procedure  used,  along  with  the  Fortran 
program,  is  given  in  Appendix  B.  The  detailed  derivation  of  the  co¬ 
state  equation,  as  needed  for  the  numerical  computation,  is  given  in 
Appendix  A.  The  final  form  of  the  costate  equation, arrived  at  in  the 
appendix, is,  however,  given  here. 


The  maximum  number  of  delays  in  the  state  will  be  2^, 


because  the  mean  incubation  period  y2  is  larger  than  the  mean  latent 
period  y-j .  Thus  each  costate  equation  requires  2y^  equations  to 
represent  it  over  the  entire  period  from  final  time  to  initial  time. 

As  derived  in  Appendix  A,  however,  each  of  these  sets  of  2y^  equations 
has  been  condensed  into  one  equation,  and  the  resulting  four  costate 


equations  are: 
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X1  "  x^  [ui(t_Tc)  +  u2(t)  +  B  y2(t)]  -  RLMD  •  y2(t)  (4.23) 


X2  =  X1  3  +  X2  u2(t)  -  2  x4  P2  *  y2(t)  "  RLMD  *  y](t) 

(4.24) 

*3  =  0  (4.25) 
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(4.27) 


Here  m(t)=o  at  final  time  t^  and  progressively  increases  by  1  each 
time  solution  proceeds  one  day  backwards  from  final  time. 


This  form  of  the  costate  equation  is  very  suitable  for  the 
numerical  solution  of  the  problem,  since  the  costate  equations  have  to 
be  solved  backwards  in  time  following  a  solution  of  the  state  equations 
in  the  forward  direction.  The  initial  values  of  the  costates  for  back¬ 
ward  integration  (actually,  values  of  costate  at  final  time)  are 
obtained  from  equations (4.21 )  and  (4.10).  These  values  are: 
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^1 ( t^)  -  0 

(4.28) 

*2^f )  ~  C-j 

(4.29) 

A3(tf)  =  C-j 

(4.30) 

^(tf)  =  1 

(4.31) 

With  the  final  values  of  costates  known,  RLMD  can  he 
easily  computed  at  the  start  of  the  backward  integration  and  is 
progressively  updated  as  the  solution  proceeds. 

4.5  PARTAN  Gradient  Algorithm 

To  obtain  the  control  function  which  minimizes  the  cost 
function,  the  negative  gradient  method  of  Kelly  [20],  coupled  with 
parallel  tangent  technique  (PARTAN)  [31],  has  been  successfully  used. 
Usefulness  of  rapid  convergence  methods  for  Min-H  problems  has  already 
been  demonstrated  by  Gottlieb  H2]  and  Lasdon  et  al  [25],  whereas  the 
steepest  descent  method  has  been  used  by  McAulay  [27]  for  systems  with 
finite  number  of  delays.  The  conjugate  gradient  method  of  Lasdon  et 
al  [25],  when  used  for  the  present  problem,  gave  good  results  but  was 
abandoned  in  favour  of  the  simpler  and  more  general  PARTAN  technique 
discussed  at  some  length  by  Pierre  [31]  and  Shah  et  al  [38].  Negative 
gradient,  satisfying  the  conditions  of  Hamiltonian  with  a  delayed 
function  of  control,  is  as  follows. 
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Let  g(u)  =  [g i ( u )  g2 ( u ) 1  be  the  gradient  vector  for  the 
two  controls,  where  g-j  is  the  gradient  of  H  with  respect  to  control  u-| 
and  g^  is  the  same  with  respect  to  control  u^. 


Then 


g(u) 


9H 

9y 


r  3  H  3H 
L3U-j  3U2J* 


for  a  no  delay  case.  As  per  the  conditions  laid  down  in  equation 
(4.22),  for  a  delay  in  control  u-|  only, 


t-+T 


Mu>  +  C’  ^  ^  1  ^  VTc  (4'32a) 


and 


3H_ 

3Un 


tf"Tc  <  t  -  tf 


(4.32  b) 


(4.33) 


Thus  we  obtain  the  gradient,  as  a  function  of  time,  as  follows: 

g-j(u)  =  [A  +  2B  u^t)]  x4  -  X1(t+Tc)  •  y-|(t+Tc),  to<t<tf-xc  (4.34) 


=  [A  +  2B  (t)]  A4, 


tr_T  <  t  <  tr- 

f  c  —  f 


(4.35) 


g^(u)  -  [C  +  2D  u 2 ( t ) U  x4  -  A-j  y-j  -  x^  y2* 


(4.36) 
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The  gradient  vector  g(u),  calculated  from  expressions 
(4.34)  to  (4.36),  is  an  implicit  funct  ion  of  time  by  virtue  of  u  being 
a  function  of  time.  Rased  on  the  algorithm  proposed  by  McAulay  [27] 
and  Lasdon  et  al  [25], for  steepest  descent  and  conjugate  gradient 
solution,  respectively,  of  such  problems,  the  following  procedure  for 
the  application  of  negative  gradient  to  our  problem  is  proposed: 

Assuming  zero  control  and  known  initial  conditions,  the 
state  equations  (4.1)  to  (4.4)  are  integrated  from  an  initial  time  t 
to  a  final  time  t^  which  satisfies  the  final  conditions  ^=0).  This 
value  of  t^  is  now  taken  as  the  terminal  time  for  subsequent  inte¬ 
grations  of  the  state  equations.  Cost  J  is  evaluated  using  equation 
(4.10)  and  the  final  values  of  the  state  already  obtained.  Thus  we 
get  the  cost  of  the  epidemic  with  no  control  applied,  and  our  aim  is 
to  reduce  this  cost  to  a  minimum  by  the  application  of  suitable  control. 
This  optimum  control  is  generated  by  the  use  of  the  following  iterative 
procedure. 

(a)  With  the  final  values  of  costates  given  by  expressions 

(4.28)  to  (4.31),  and  using  the  stored  values  of  the  states,  just 
calculated  (for  the  entire  period  of  time),  the  costate  equations 
(4.23)  to  (4.26)  are  solved  by  backward  integration  from  final  time 
tf  to  initial  time  tQ.  Thus  we  have  the  values,  necessary  for  the 
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evaluation  of  gradient  vector  g(u).  If  the  control  vector  for  the 
i  th  iteration  is  represented  by  u1  and  the  corresponding  gradient 
by  g(un ) ,  let  us  assume 

s1  =  -  gfu1).  (4.37) 

s1 ,  in  its  present  form,  represents  the  negative  gradient  for  the  i  th 
iteration,  and  has  two  components,  s-j  and  s^»  corresponding  to  two 
controls,  u-j  and  u^.  Reckoning  of  the  iteration  number  i  may  start 
either  from  0  or  1 ,  depending  upon  the  personal  preferences.  In  the 
present  case,  let  us  assume  i  to  start  from  0. 


(b)  We  now  proceed  to  search  for  a  multiplying  factor  a 

consisting  of  two  components  ct-|  and  c^*  such  that  a  value  of  increment 
control  (in  the  vector  sense)  may  be  obtained  to  minimize  cost  J.  Thus 
choose , 


i  .  ...  .  /  i  ,  i  i  \ 

a  =  a  to  minimize  J(u  +  ct  s  ). 


(4.38) 


Then  the  new  control  vector  will  be. 


•  .  T  •  •  • 

1+1  1  ,  11 

U  =  U  +  a  S  . 


(4.39) 

The  two  components  of  a  are  evaluated  by  conducting  a  one  dimensional 
search,  first  for  one  control  and  then  for  the  other,  to  minimize  J  in 
each  case. 


(c)  Using  the  new  control  vector  obtained  from  equation 
(4.39),  we  proceed  to  solve  the  state  equations,  anew,  and  calculate 
the  new  value  of  cost  J.  If  this  new  value  of  J  is  judged  to  be  the 


minimum  value  the  solution  can  stop  here,  otherwise  we  proceed  to  the 
next  step. 

(d)  Now  we  either  calculate  the  new  costates  and  new 

gradient  using  the  procedure  outlined  in  step  (a)  or  find  the 
acceleration  step,  and  go  to  step  (b).  The  decision  as  to  whether 
this  step  is  to  be  a  steepest  descent  or  an  acceleration  step  is  made 
on  the  basis  of  PARTAN  procedure  discussed,  at  some  lenqth,  by  Pierre 
[31]  and  Shah  et  al .  [38],  and  outlined  below: 

Controls  u  and  u  are  obtained,  respectively,  from  u° 
and  u1  by  using  the  steepest  descent.  An  acceleration  step,  which  is 

r\ 

the  vector  difference  of  controls  u  and  u°  is  now  obtained  and 

treating  it  as  negative  gradient,  control  u  is  obtained.  Controls 

U4,  u6,  u8  ...  are  now  obtained  from  controls  u8,  u8,  u7  ..., 

respectively,  by  the  use  of  steepest  descent  approach;  whereas  u  ,  u  ... 

etc.  are  obtained  from  the  corresponding  pairs  of  controls  u4,  u^  and 
6  3 

u  ,  u  ...  etc.,  by  the  use  of  acceleration  step. 

The  process  is  iterated  until  the  cost  function  J  converges 
to  a  minimum  value.  The  numerical  procedure  based  on  this  algorithm 
is  given  in  Appendix  B.  It  may  be  pointed  out  here  that,  although  the 
gradient  method  assumes  an  arbitrary  6u ,  there  is  a  constraint  in  the 
present  case.  Since  a  negative  control  is  unthinkable  in  practice 
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(vaccinations  cannot  be  undone!),  while  searching  the  values  of  a 
the  values  of  trial  control  which  would  otherwise  be  negative  were 
replaced  by  zero.  The  final  values  of  cost  obtained  are  therefore, 
minimum  in  the  practical,  constrained  case.  The  results  of  the 
numerical  solution,  presented  in  the  next  section,  prove  the  effect¬ 
iveness  of  the  method. 

4.6  Choice  of  cost  constants 

While  formulating  the  optimum  control  problem  in  section 
4.2,  the  cost  constants  for  immunization  controls  were  represented 
by  A,  B,  C  and  D,  and  that  for  the  cost  per  reported  case  by  C-j . 
Factors  influencing  the  choice  of  numerical  values  for  these  constants 
are  discussed  in  this  section;  a  choice  is  also  made  for  a  represent¬ 
ative  set  of  values  to  be  used  for  the  numerical  solution  of  the 
problem. 

Since  we  are  using  the  normalized  variables  in  the  model, 
the  state  variables  y-j  ,  y^»  y^  and  y^  represent  per  unit  values.  The 
actual  values  are,  therefore,  obtained  by  multiplying  these  values  by 
the  reference  value,  which  is,  in  this  case,  the  number  of  individuals 
forming  the  population  of  the  community.  For  y^  to  have  this 
dimension,  constants  A  and  C  should,  respectively,  represent  the  cost 
of  one  inoculation  of  active  and  passive  immunization  agents. 
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Similarly,  the  cost  function  J  will  also  be  dimensionally  consistant 

if  represents  the  cost  of  one  individual  falling  sick  with  the 
di sease. 

Although  the  cost  of  immunization  per  inoculation  is  fairly 
constant  when  such  inoculations  are  given  at  a  slow  rate,  yet  it  is 
far  from  realistic  to  consider  this  cost  to  be  constant  in  epidemic 
situations.  In  such  situations  the  supply  of  vaccine  or  antiserum, 
as  the  case  may  be,  has  to  be  supplemented  by  creating  additional 
manufacturing  capacity,  and  crash  immunization  programs  have  to  be 
initiated,  thus  causing  the  cost  per  injection  to  rise.  Although  the 
full  effect  of  this  dislocation  of  normal  services,  on  the  per  unit 
cost  of  immunization,  is  not  yet  known  precisely,  yet  it  is  certain 
that  the  cost  per  unit  of  immunization  is  a  non-linear  function  of 
the  rate  of  immunization.  The  simplest  way  of  representing  a  non¬ 
linear  quantity  is  to  represent  it  as  the  first  two  terms  of  a  Taylor 
series  expansion.  This  fact  is  utilized  here  and,  as  such,  the  cost 
of  active  immunization  is  assumed  to  be  A  +  Bu-j  per  injection  and  that 
of  passive  immunization  to  be  C  +  Du^  per  antiserum  injection.  When 
proceeding  with  the  optimization  process,  constants  B  and  D  penalize 
high  control  effort,  hence  prevent  the  optimum  control  rates  becoming 
unrealistically  high. 
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Accurate  values  of  constants  A,  B,  C  and  D  can  only  be 
determined  after  a  sufficient  experience  with  the  model  and  a  full 
knowledge  of  the  effect  of  immunization  rates  on  the  unit  cost  of 
immunization.  However  it  can  be  generally  assumed  that  the  cost  of 
one  injection  of  antiserum  is  relatively  higher  than  that  of  the 
vaccine  for  the  same  disease,  and  that  B/A  and  D/C  are  relatively 
large  values.  The  following  representative  values  have  been  assigned 
to  these  constants,  for  the  numerical  solution  of  the  problem. 

A  =  $  2.oo,  B  =  $  200. oo 

C  =  $  5.oo,  D  =  $  500. oo. 

This  choice  of  relative  values  for  linear  (A  and  C)  and  quadrature 

(B  and  D)  cost  constants  is  quite  consistant  with  the  economics  of 

vaccine  and  antisera  production,  in  practice.  Restricting  our  attention 

to  active  control  only,  we  can  say  that  for  very  small  values  of  u-j 

the  cost  per  inoculation  is  approximately  equal  to  A,  whereas  it  is, 

almost,  independent  of  A  for  very  large  values  of  u^  A  value  of  u-j 

A  1 

for  which  A  and  B  have  equal  weightage  is  ^  per  day  or  3.65  per 

year.  This  corresponds  to  a  production  capacity  which  can  be  geared 
to  produce  3.65  vaccines  per  person  per  year.  If  a  control  at  a  higher 
rate  is  required,  additional  capacity  has  to  be  created  at  a  high 
fixed  cost,  thus  resulting  in  a  relatively  high  per  vaccine  cost, 
which  is  dependent  on  B  only.  The  same  is  also  true  for  passive 


-  85  - 


control . 


Cost  per  reported  case  (C-|)  depends  on  the  disease  in 
question.  It  will  be  higher  for  a  disease  which  causes  the  victim  to 
be  away  from  work  longer.  Moreover,  it  will  be  higher  for  those 
diseases  which  are  costlier  to  cure  and  those  which  have  higher 
fatality  rates.  Three  representative  values  of  C-j  considered  in  this 
thesis  are  $  50. oo,  $  100. oo  and  $  500. oo. 

4.7  Results  of  optimization 

The  optimization  procedure  discussed  in  this  chapter  has 
been  applied  to  the  epidemic  problem  discussed  in  Chapter  III. 

Results  of  the  numerical  solution  of  the  problem  are  presented  in  this 
section.  For  the  sake  of  clarity  the  parameters  of  the  control 
problem  are  restated  below: 


Epidemic  Parameters 

Initial  susceptibles  (y-j(t0)) 
initial  infectives  (y2 (t  ) ) 

y2(t),  t  <  to 
effective  contact  rate  (b) 
factor  k(t) 

mean  value  of  the  latent  period 
mean  value  of  the  incubation  period 


=  0.5, 

=  10"4, 

=  0.0, 

=  1.5, 

=  0.9, 

=  7  days, 

=  1 4  days , 
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standard  deviation  of  the  latent  period 

=  1 .3, 

standard  deviation  of  the  incubation  period 

=  2.3. 

Control  Parameters 

Per  unit  cost  of  active  control  (A) 

=  2.0, 

square  law  cost  of  active  control  (B) 

=  200.0, 

per  unit  cost  of  passive  control  (C) 

=  5.0, 

square  law  cost  of  passive  control  (D) 

=  500.0, 

penalty  function  (p^) 

=  1 06 , 

cost  per  reported  case 

=  $  50.0, 

$  100.0  and 
$  500.0. 

The  choice  of  these  parameters  has  already  been  discussed 
in  the  various  sections  of  this  thesis.  It  may,  however,  be 
emphasized  once  more  that  these  parameters  do  not  represent  any 
particular  disease,  but  are  typical  values  used  to  demonstrate  the 
usefulness  of  the  procedure  developed  in  this  thesis.  Actual  values, 
for  a  particular  case,  can  be  used  to  get  the  results  appropriate  to 
that  application. 

In  addition  to  the  case  where  the  active  control  is 
assumed  to  be  immediately  effective  after  it  is  administered,  two  cases 
of  delay  in  active  control  have  also  been  considered.  The  delays 
considered  are  5  days  and  16  days,  the  latter  being  larger  than  the 
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mean  incubation  period.  Thus,  in  all,  nine  cases  have  been  considered. 
The  results  of  numerical  solution  of  these  nine  cases  are  summarized 
below: 

Results 

The  solution  of  the  model,  in  all  the  nine  cases,  and  when 
no  control  is  used,  is  the  same  as  that  shown  by  the  no-control  curve 
in  figures  3.1  to  3.7.  The  terminal  time  is  found  to  be  145  days  in 
each  case  and  the  final  extent  of  the  epidemic,  without  control,  is 
0.4457  reported  cases  (per  unit  value).  Therefore  t^  used  for  the 
optimization  problem  is  145  days. 

Figure  4.1  shows  the  convergence  of  the  optimization 
process  to  an  optimum  cost,  for  a  typical  case,  whereas  the  results 
relating  to  all  the  nine  cases  considered  are  summarized  in  table  4.1. 
This  table  shows  the  optimum,  in  each  case,  calculated  against  the 
cost  without  control.  Optimum  costs  expressed  as  percentage  of  the 
initial  cost,  and  the  percentage  saving  in  these  costs  due  to  the  use 
of  optimum  control  strategy,  are  shown  in  the  next  two  columns. 

Column  6  shows  the  final  extent  of  the  epidemic,  in  each  controlled 
case,  as  a  percentage  of  that  without  control. 

The  optimum  active  and  passive  controls,  determined  by  the 
optimization  procedure,  are  plotted  in  figures  4.2,  4.3  and  4.4. 


COST  IN  DOLLARS 
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ITERATION  NUMBER 


FIG.  4-1 


CONVERGENCE  OF  THE  COST  TO  OPTIMUM  VALUE. 
COST  PER  REPORTED  CASE  -  $  500.00,  AND  NO 
DELAY  IN  ACTIVE  CONTROL. 
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table  4-1 
COST  CONSTANTS 

A  =  2.0  _  B  =  200.0  C  =  5.0  D  =  500.0 


COST  PER 
REPORTED 
CASE  IN  $ 

COST 

WITHOUT 

CONTROL  IN  $ 

OPTIMUM 

COST  IN  $ 

OPTIMUM 

COST  AS  % 

OF  COST 

WITHOUT 

CONTROL 

%  SAVING 

IN  COST,  BY 

OPTIMIZATION 

REPORTED 

CASES  WITH 

CONTROL  AS 

%  OF  THOSE 
WITHOUT 

CONTROL 

DELAY  IN 

ACTIVE 

CONTROL 

50.0 

$  22,290 

11.518 

51.7 

48.3 

18.05 

NO  DELAY 

12.109 

54.3 

45.7 

19.55 

5  DAYS 

13.703 

61.5 

38.5 

28.40 

16  DAYS 

100.0 

$44,581 

13.467 

30.2 

69.8 

7.15 

NO  DELAY 

14.867 

33.4 

66.6 

8.14 

5  DAYS 

17.849 

40.05 

59.95 

12.35 

16  DAYS 

500.0 

$  222.904 

21.489 

9.61 

90.39 

1.95 

NO  DELAY 

24.121 

10.8 

89.2 

2.34 

5  DAYS 

29.998 

13.4 

86.6 

3.47 

16  DAYS 
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Figure  4.2  shows  the  optimum  controls  for  three  different  case  report 
costs,  without  any  delay  in  the  effect  of  active  control.  Figures 
4.3  and  4.4  show,  respectively ,  the  corresponding  controls  for  the 
cases  of  6  days  delay  and  16  days  delay  in  active  control. 

4.8  Discussion  of  the  results 

The  results  presented  in  the  last  section  clearly 
establish  the  usefulness  of  the  optimization  procedure  developed  in 
this  thesis.  It  is  apparent  from  figure  4.1  that  the  cost  function 
converges  to  a  minimum  value  (which  is  substantially  lower  than  the 
initial  cost)  rapidly.  A  similar  decrease  in  cost  is  also  observed 
in  the  other  cases  tabulated  in  table  4.1.  We  find  that  the  saving  in 
the  total  cost  ranges  from  38.5%  to  90.39%.  The  percent  saving  in 
cost  is  highest  in  the  case  of  the  most  severe  disease,  with  cost  per 
reported  case  of  $  500. oo,  because  in  this  case  it  is  economical  to 
use  more  control  effort  to  prevent  the  disease.  This  fact  is  also 
demonstrated  by  the  observation  that  the  actual  disease  cases  have 
been  reduced  to  less  than  4%,  compared  to  over  28%  in  the  first  case 
(where  C-|  =  $  50. oo).  Another  important  observation  is  that  the 
saving  in  cost,  by  optimization,  is  reduced  when  delay  in  the  active 
control  is  introduced.  These  results  are  in  line  with  the  intuitive 
thinking  that  more  control  should  be  applied  for  a  disease  which 
causes  higher  losses  and  that  delay  in  control  reduces  its  effectiveness. 


PASSIVE  CONTROL  u,  ACTIVE  CONTROL 
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TIME  (days) 

FIG.  4-2  OPTIMUM  CONTROL  PLOTTED  AGAINST  TIME. 

NO  DELAY  CASE. 


PASSIVE  CONTROL  u9  ACTIVE  CONTROL 


-92- 


FIG.  4-3  (a)  ACTIVE  CONTROL  WITH 


FIG.  4-3  OPTIMUM  CONTROL  PLOTTED  AGAINST  TIME. 
DELAY  OF  5  DAYS  IN  ACTIVE  CONTROL. 


PASSIVE  CONTROL  u,  ACTIVE  CONTROL 
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FIG.  4-4  (a)  ACTIVE  CONTROL  WITH 


TIME  (days) 

FIG.  4-4  OPTIMUM  CONTROL  PLOTTED  AGAINST  TIME. 
DELAY  OF  16  DAYS  IN  ACTIVE  CONTROL. 
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Turning  to  the  plots  of  optimum  control,  we  find  that 
larger  amounts  of  control  are  used  as  the  cost  per  case  increases 
from  50. oo  to  500. oo.  Moreover,  we  find  that  more  active  control  has 
been  used  than  passive  control,  which  is,  perhaps,  due  to  the  relative 
cheapness  of  the  active  control.  We  also  find,  as  one  would  expect, 
that  with  delay  introduced  in  the  active  control,  the  use  of  passive 
control  is  increased. 

Thus  we  find  that  the  optimization  procedure  developed 
in  this  thesis  gives  encouraging  results  that  are  consistant  with 
i ntuit ion. 

4.9  Optimum  control  procedure  for  case  2 

We  have  so  far  restricted  our  attention  to  case  1  only. 

This  was  done  to  keep  the  discussion  unambiguous.  The  method  is, 
however,  equally  applicable  to  case  2.  Using  the  procedure  developed 
in  section  4.4  and  detailed  in  Appendix  A,  the  corresponding  costate 
equations  for  case  2  are  as  follows: 

A.|  =  X1  {u-|(t-Tc)  +  6  y2(t)>  “  RLMD  *  >'2^)  (4.37) 

X2  =  x1  6  y-j  (t)  -  2x4  p2  «(y2)  •  y2(t)  -  RLMD  •  y-,  (t)  (4.38) 

x3  =  0  (4.39) 

x4  =  0, 


(4.40) 
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where 

m(t) 

RLMD  =  ^  WLP(n)  •  X2(t+nx)  -  WIP(n)  r*2(t+riT)  ^  +  u2(t+nx)} 

n=o 

-  X^(t+nx)  -  x^(t+m)  {C  U2(t+nx)  +  D  U2^(t+nr)}],  (4.41) 

m(t)  is  the  same  integer  variable  as  in  equation  (4.27). 

The  corresponding  gradient  equations  are: 

g-j(u)  =  [A  +  2B  u](t)]x4  -  X1(t+xc)  •  y-|(t+xc),  tQ  <  t  <_tf-xc 

=  [A  +  2B  u1(t)]x4,  tf-xc  <  t  <  tf  (4.42) 

2y  2 

g2(u)  =  [x4(C+2  U2(t)-D)-X2(t)]  WIP(n).y1(t-nx)-y2(t-nx).(4.43) 

n=o 

When  the  same  values  of  constants  were  used  in  this  case, 
substantial  reductions  of  cost  were  observed.  The  final  cost  was, 
however,  found  to  be  relatively  less  sensitive  to  passive  control  which 
is  due  to  the  fact  that  the  nature  of  U2  here  is  different  than  that  in 
the  previous  case.  U2,  here,  is  the  number  of  passive  immunizations 
given  per  reported  case,  and  not  the  actual  number  of  immunizations. 

The  choice  of  C  and  D  should,  therefore,  be  revised  to  suite  this 
case.  A  reference  to  equation  (4.8)  shows  that  the  nature  of  cost 
constant  C  is  substantially  same  as  that  in  case  1 ,  and  it  is  the 
cost  of  one  antiserum  injection.  But,  the  choice  of  constant  D  is 
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more  complicated.  D,  in  the  present  case,  not  only  represents  the 
increase  in  the  cost  per  inoculation  due  to  the  increase  in  inocul¬ 
ation  rate,  but  also  represents  the  increase  in  cost  due  to  the 
added  cost  of  more  intensive  search  for  the  contacts  of  the  reported 
case.  Effect  of  the  two  factors  on  the  choice  of  D  is  not  fully 
known  and  more  study  is  needed  for  a  better  estimation  of  D.  It  is 
however  certain,  that  with  a  better  choice  of  these  constants,  useful 
results  can  be  obtained  with  this  model  also. 

4. 10  Conclusions 

We  thus  conclude  that  the  optimization  procedure  developed 
in  this  chapter  is  quite  effective  in  determining  the  optimum  control 
of  an  epidemic.  Although  the  continuously  changing  rate  of  control, 
as  observed  in  the  results  presented  in  section  4.7,  would  be 
difficult  to  implement  in  practice,  yet  it  serves  as  a  guide  for  the 
control  strategy  to  be  used  in  practice.  If  the  control  strategy,  in 
practice,  has  to  be  varied  due  to  some  other  considerations,  the 
quantitative  effects  of  this  deviation  can  be  easily  assessed  using 
this  model.  This  will  help  provide,  at  least,  a  quasi-optimal  control 
strategy  in  that  case. 
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CHAPTER  V 

SUMMARY  AW  CONCLUSIONS 

5.  I  Summary 

A  mathematical  model,  which  can  be  used  for  the  prediction 
and  optimum  control  of  all  contagious  diseases  that  are  transmitted 
through  personal  contact  between  infectives  and  susceptibl es ,  has 
been  developed.  The  new  deterministic  model  takes  into  full  consider¬ 
ation,  for  the  first  time,  the  latent  and  infections  periods  of  the 
disease  and  their  statistical  variations.  The  model,  though  only 
applicable  to  closed  populations,  yet  gives  results  which  are  nearer 
to  disease  spread  observed  in  practice  than  those  obtained  from 
earlier  deterministic  models.  A  hypothesis  about  the  possible 
relation  of  meteorological  changes  to  disease  outbreaks,  based  on  the 
results  obtained  from  the  new  model,  has  also  been  given. 

The  sensitivity  of  the  model  to  active  and  passive  immu¬ 
nization  controls  has  been  studied,  and  a  procedure  for  charting  an 
optimum  control  strategy,  by  the  application  of  Pontryagin's  Minimum 
Principle,  has  been  developed.  For  the  first  time,  dynamic 
optimization  techniques  have  been  successfully  used  for  finding  the 
optimum  control  strategy  for  this  class  of  problems.  In  addition  to 
being  a  practical  application  of  the  existing  control  theory  to  the 
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control  of  epidemics,  the  system  discussed  also  forms  an  important 
example  of  the  application  of  this  theory  to  a  non-linear  case  with 
a  finite  number  of  delays,  both  in  state  and  control.  The  results 
obtained  show  the  sensitivity  of  the  control  to  delays  in  active 
control  and  to  the  severity  of  the  disease  (represented  by  cost  per 
reported  case). 

5. 2  Conclusions 

From  the  analysis  given  in  this  thesis,  we  can  now  conclude 
that  a  contagious  disease  can  be  generally  represented,  mathematically, 
by  three  important  parameters,  viz.  effective  contact  rate,  latent 
period,  and  infective  (or  incubation)  period.  Effective  contact  rate 
(f),  the  first  of  the  three  parameters,  is  as  much  dependent  on  the 
social,  hygienic  and  environmental  conditions  of  the  population  as  on 
the  infectiousness  of  the  disease,  and  is  proportional  to  the 
infectiousness  of  the  disease,  other  conditions  being  equal.  Another 
important  conclusion,  from  the  solution  of  the  model,  is  that  for  a 
disease  to  become  epidemic,  both  the  proportion  of  initial  susceptibles 
in  the  population  and  6  should  be  greater  than  their  respective 
threshold  values. 

In  the  area  of  control  application,  we  conclude  that  the 
available  control  must  be  applied  as  early  as  possible  and  at  the 
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highest  rate.  Further,  if  the  relative  costs  of  administering  each 
control  are  known, it  is  possible  to  find  the  best  control  strategy 
by  the  use  of  optimum  control  theory. 

5.3  Recommendations  for  further  research 

This  thesis  presents  many  possibilities  for  future  research. 
The  first  necessity,  for  getting  practically  applicable  results,  is 
the  accurate  assessment  of  parameters.  Statistical  estimation  of 
these  parameters  can  be  done  more  accurately  if  the  data  of  the  right 
kind  are  available.  It  will,  therefore,  be  advisible  to  record  the 
disease  data  with  the  new  models  in  view  and  develop  methods  for 
estimation  of  the  parameters.  Laboratory  research  about  infectious¬ 
ness  of  the  diseases,  their  latent  periods,  and  their  incubation 
periods  can  supplement  statistical  research.  More  research  is  also 
needed  to  evaluate  the  cost  constants  for  the  cost  function,  more 
realistically.  Introduction  of  other  controls,  not  represented  in 
the  present  study,  will  help  make  the  model  more  useful. 

The  present  model  was  restricted  to  closed  populations. 
Extension  of  the  optimization  approach,  used  here,  to  models  modified 
to  include  parameters  relating  to  births,  deaths  and  population 
migration,  would  make  an  interesting  study.  Moreover  the  present 
study  was  restricted  to  short  term  control  of  one  outbreak;  its 


-100- 


extension  to  endemic  diseases  on  the  long-term  basis  is  another 
challenging  area. 

All  individuals  have  been  considered  to  have  similar 
response  to  the  disease.  The  effect  of  age  groups,  sex  groups  or  any 
other  population  classification  remains  to  be  studied. 

A  very  important  factor,  not  considered  here  but  already 
voiced  in  the  literature  [2],  [4],  is  the  question  of  geographical 
spread.  In  this  study  the  population  has  been  lumped  as  a  single 
homogeneous  group,  which  is  quite  far  from  reality  when  large  areas 
are  considered.  Application  of  optimization  techniques  to  models  for 
geographical  spread  may  require  the  use  of  distributed  parameters, 
and  is ,thus ,1 i kely  to  prove  a  challenging  but  useful  field  of  study. 

Only  a  complete  analytical  solution  of  the  control  problem 
can  really  establish  the  limitations  and  full  potential  of  the 
optimization  techniques  used;  obtaining  these  solutions  is  another 
area  where  considerable  research  is  yet  needed. 

If  this  thesis  helps  generate  some  interest  in  the  above 
areas  of  research,  it  will  have  served  a  purpose. 
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APPENDIX  A 

This  appendix  details  the  derivation  of  costate  equations 
and  the  assumptions  made  to  arrive  at  the  results. 

With  the  assumption  and  definitions  of  Chapter  IV,  the 
state  equations  of  the  model  can  be  written  as: 


y-i  =  -  e  y-i  y2  -  u1  (t-xc)  •  y-,  -  u2  (t)  •  y1 

2yl  r 

^  exp  {•  2  ( v1 1 —  ht )  } 

/  CT 

1  n=o  L 


( A- 1 ) 


y2  = 


_  KB  J 

/  2tt  ct 


2o 


y-|(t-nT)  •  y2(t-nT) 


KB  1 


2y, 


—  V  exp  { - 2  ( U 2~ tir )  ) 

"2  ^0  L  2o2 


y-,  (t-nr)  •  y2(t-nx) 


-  u2(t)  •  y2 


y^  = 


_  Kb  1 


2  U. 


/  2tt 


'2  n=o 


r— "  J  ^ 

\  exp  { - 2  (y2-ni)  } 

2a„ 


(A-2) 


y-j  (t-nx)  •  y2(t-nx)  (A-3) 


y4  =  «(y2)  +  A  U-,  (t)  +  B  u2(t)  +  [C  u2(t)  +  D  u2(t)]  ,  (A-4) 


where 


«  (y2)  =  o 


y2  >  o 


<s  (y2)  =1  y2  <  o, 


and  p2  is  the  penalty  multiplier. 
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The  constant  multipliers  and  exponential  factors  in  the 
above  expressions,  the  latter  representing  the  distribution  of  latent 
and  incubation  periods,  can  be  calculated  and  stored  once  for  all, 
for  all  values  of  n. 


Substituting  the  stored  vectors  WLP(n)  and  WIP(n),  defined 


as 


WLP(n) 

 Ke 

1_ 

exp  { - — 9 

(u-.-nx)^} 

/  2tt 

°1 

2a/ 

1 

WIP(n) 

_  K3 

/  2xr 

1_ 

a2 

exp  { — — ? 
2  a2 

(u2~nx)^} 

the  state  equations  become: 

Yl  =  -  3  y-|  y2  -  ^(t-x^  *  yl  ‘  u2^  *  yl 

2y i  2^2 

y2  =  Y_  WLp(n)  *  y]  (t-m)  •  y2(t-nx)-  WIP(n)  x 
n=o  n=o 


y-j  (t-nT)  •  y2(t-nx) 


2y. 


y3  =  WIP(n)  •  y-j(t-nx)  •  y2(t-nT) 


( A- 5 ) 


-  u2(t)  •  y2  (A-6) 


(A-7 ) 


n=o 


y4  =  p2  y2  S(y2)  +  A  U](t)  +  B  u/t)  +  C  u2(t)  +  D  u|(t)  (A-8) 


This  version  of  the  state  equations  is  quite  suitable  for 
numerical  analysis  as  they  can  be  easily  converted  to  the  differential 
difference  form.  The  equations  can  be  rewritten  in  a  concise  form  as 
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below: 

y-|  =  -  b  y-,  y2  -  u1(t-T(.)  •  y,  -  u2(t)  .  y1  (A-9) 

y2  =  ACTV  -  RMVD  -  u2(t)  •  y2  (A-10) 

y3  =  RMVD  (A-ll) 

y4  =  p2y2  «(y2)  +  A  u^t)  +  B  u^(t)  +  C  u2(t)  +  D  u|(t)  (A-12) 


where 


2y- 


Y  WLP(n) 


n=o 


y-j  (t-nT) 


y2(t-nT)  =  ACTV 


(A-l 3a) 


is  the  rate  at  which  latent  cases  become  active  in  spreading  the 
disease,  and 


n=o 


(A-l 3b) 


is  the  rate  of  reported  cases  per  unit  time,  at  time  1 t ' .  Since  the 
unit  of  time  used  is  one  day,  the  above  rates  can,  respectively,  be 
called  the  per  day  active  and  removal  rates. 
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COSTATE  EQUATIONS 


For  a  complete  representation  of  the  process  of  disease 
spread,  our  model  needs  the  past  history  of  state  variables  for  twice 
the  incubation  period.  Since  the  values  of  state  variables  have  been 
stored  at  intervals  of  one  day,  a  total  of  delayed  values  of  each 
variable  may  figure  in  the  state  equations.  As  per  the  conditions 
for  formulation  of  costate  equations,  arrived  at  by  McAulay  [27]  and 
others  [9],  [22],  [32]  and  discussed  in  Chapter  4,  a  total  of  2y^ 
expressions  are  needed  to  completely  represent  each  costate  between 
the  initial  and  final  time.  Working  backwards  in  time,  2^-1  of 
these  expressions  represent  the  costate  in  question,  one  for  each  time 
slot  of  one  day  duration  between  final  time  t^  and  t^-2y2.  Another 
one  is  needed  for  the  remaining  time  up  to  the  initial  time. 


The  costate  equations  for  the  last  interval,  or  the  first 
one  starting  backwards  from  final  time  t^,  can  be  written  from  the 
general  formula 


x 


X 


X 


1 

2 

3 


ML  = 

avi 

ML  = 
9y2 


af 1  a  f ^  af3  9f4 

3y7  "  X2  "  A3  3^7  '  ^4  3^ 
af 1  9  f  ^  9f  3  9f4 

X1  ay^  ‘  x2  ay^  '  A3  ay^  '  x4 

af1  af2  9f  3  sf4 

X1  3^'  x2  '  x3  ay^  '  x4  8^ 


(A-14) 
(A— 1 5) 

(A-16) 
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•  _  9H 

A4  '  '  W, 


9fn 


=  -  X 


lay/ 


-  x 


9f , 
2  97 


-  X 


9f. 
3  97 


-  X 


3f^ 

4  9y, 


(A-l  7) 


When  applied  to  the  present  problem,  it  is  apparent  that 
the  right  hand  sides  of  equations  (A-l 6 )  and  (A-l 7 )  are  zero,  so  only 
the  first  two  equations  remain  to  be  considered.  Concentrating  on 
equation  (A-l 4 )  we  find  that  this  equation,  between  time  t^-i  and 
t.f-2 t,  will  be 


•  =  9H_  9H  | 

1  "  3y1  "  9y-|  (t-x)  ]t=t+T. 


(A-l 8) 


Equation  ( A- 18)  is  the  same  as  equation  (A-l 4)  with  one  term  added  to 
its  right  hand  side.  This  additional  term  is  the  partial  derivative 
of  Hamiltonian  with  respect  to  the  delayed  variable  y-|(t-x)  and  the 
resulting  term  advanced  by  t  in  time.  Thus  an  equation  for  the  time 
between  tf-kx  and  tf-(k+l)x,  will  have  k  similar  terms  added  to  the 
right  hand  side  of  (A-l 4) -  Hence  the  corresponding  costate  equation 
for  time  t^-kx  <  t  <  t^-(k+l)x  is: 

.  _  9H  9H _  ,  _  _  9H  | 

X1  "  3y-|  "  9y i  (t-x)  't=t+x*’*'  9.y-|  ( t-kx )  't=t+kx  (A-l 9) 


Rewri ti ng 


9H 

sy-i 


k 

L  8y-|  ( t-rvrT  U't+nr*  k  =  1-2’3'  "•  Zv2 ' 
n=l  1  (A-20) 


Similarly 


k 


9H 

9y2 


n=l 


9H 

3y2(t-nT) 


U=t+nt’  k  =  1  ,2,3,  2u2' 

(A-21) 
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Substituting  the  expression  for  Hamiltonian  and  taking 
necessary  partial  derivatives,  the  costate  eguations  for  time 
t_p-kx  <  t  <  t^-(k+l)x  become: 

k 

X1  =  6  y2  X1  +  X1  ul(t'Tc'  +  X1  “2^)  '  y2(t)  Y  X3(t+nT)  •  WIP(n) 


n=o 


-  y2(t)  ^2^t+nx^  '  tWLP(n)  -  WIP(n)] 


n=o 


where  k  =  0,  1 ,  2  _  2y, 


(A-22) 


A2  =  6  yl  A1  +  a2  u2^  "  2x4  P2y2  6^y2^  "  yl^'  Y  X3^t+nT^  *WIP(n) 


n=o 


-  y^t)  y  X2(t+nx)  •  [WLP(n)  -  WIP(n)] 


n=o 


where  k  =  0 ,  1 ,  2  ....  2y^ 


(A-23) 


X3  ^ 


X4  =  0 


(A-24) 

(A-25) 


Eguations  (A-22)  and  (A-23)  simplify  to 


X-,  =  x1[u1(t-xc)  +  u2(t)  +  6  y2(t)l  -  y2(t)  Y  Mt+nx) 


WLP(n) 


n=o 


+  y  (t)  V  WIP(n)[x2(t+nT)  -  x3(t+m)] 


n=o 


2y . 


where  k  =  0,  1  ,  2,  3 


•  •  •  • 


(A-26) 
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2  "  A1  3  yl(t)  +  a2  u2  "  2x4P2y2  6^y2)  '  yl(t)  A2(t+nx)  ‘  WLP(n) 
k 

+  yl^  ^  WIP(n) [x2(t+nx)  -  A3(t+nt)] 


n=o 


n=o 


where  k  =  0,1 ,2,3  _  2y< 


(A-27 ) 


Defining  a  variable  RLMD  as: 


k 

RLMD  =  ^  WLP(n)  •  X2(t+nx)  -  WIP(n)  •  [x^ ( t+nx )  -  x3(t+nt)] 
n=o 


where  k  =  0,1 ,2,3  ....  2y^ 


(A-28) 


The  costate  equations  reduce  to 

^  =  X-|[u-|(t-xc)  +  u2(t)  +  6  y2(t)]  -  RLMD  •  y2(t)  (A-29) 

X2  =  x-|  6  y-| ( t)  +  x2  u2(t)  -  2  x4  P2  6(y2)  •  y2(t)  -  RLMD  •  y-, ( t ) 

(A-30) 

X3  =  0  (A-31 ) 

X4  =  0  (A-32) 

Costate  equations  (A-29)  to  (A-32)  are  apparently  inde¬ 
pendent  of  time  slot  being  considered  but  it  is  not  so  in  reality.  The 
newly  defined  variable  RLMD,  which  figures  prominently  in  the  costate 
equations,  is  completely  dependent  on  the  time  slot,  and  has  to  be  up- 
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dated  at  every  step  as  the  backward  integration  of  costate  equations 
proceeds.  A  look  at  the  expression  (A-28)  reveals  that  when  k  =  0, 
the  summation  is  over  only  one  value  of  n,  and  the  value  of  RLMD 
depends,  only,  on  the  known  final  values  of  state,  costate  and  control 
variables.  The  costate  differential  difference  equations  can  now  be 
solved  for  one  step  in  time;  a  day  in  the  present  case.  Expression 
(A-28)  is  processed  again;  now  with  k  =  1.  The  value  of  RLMD  now 
depends  on  the  values  of  state,  costate  and  control  variables  between 
periods  t^  and  t^-x.  These  values  are  known  by  this  time.  The  inte¬ 
gration  proceeds  in  this  manner  raising  the  value  of  k  by  one,  each 
step  of  time,  till  it  reaches  2^.  After  this  moment  k  remains  fixed 
at  this  value  but  RLMD  is  still  updated  every  step  of  time  up  to  the 
end  of  integration  procedure  at  t=tg.  Thus  we  are  able  to  evaluate  the 
costates  for  the  entire  period  of  time. 
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APPENDIX  B 

This  appendix  outlines  the  numerical  procedure  used  for 
the  solution  of  the  optimum  control  problem  discussed  in  this  thesis. 
The  Fortran  IV  program,  suitable  for  the  model  IBM  360/67  digital 
computer,  and  successfully  applied  to  the  determination  of  an  optimum 
control  strategy,  is  given  at  the  end  of  this  appendix. 

For  the  purpose  of  numerical  simulation,  the  discretized 
version  of  the  model  given  in  Chapter  IV  and  represented  by  equations 
(4.1)  to  (4.4)  for  case  1  and  by  equations  (4.5)  to  (4.8)  for  case  2 
has  been  used.  A  unit  of  time  of  one  day  has  been  used  to  make  the 
report-rate  data,  generated  by  the  simulated  model,  comparable  to  the 
field  data.  This  unit  of  time  may,  however,  be  replaced  by  a  week  or 
a  month  in  the  case  of  slow  spreading  diseases  like  tuberculosis. 

At  the  beginning  of  any  given  day,  the  values  of  variables 
a ( t )  and  r(t)  (rates  of  newly  active  cases  and  reported  cases)  are 
calculated  using  the  stored  values  of  the  history  of  the  states  and 
the  corresponding  normally  distributed  weighting  multipliers.  These 
rates  are  now  assumed  constant  for  the  day  and  the  state  equations  are 
solved  for  one  day  using  Runga  Kutta  method.  Final  values  of  the 
states  at  the  end  of  one  day  form  the  initial  values  at  the  beginning 
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of  the  next  day.  This  process  is  carried  on  till  the  stopping 
criterion,  which  is  either  reaching  the  assigned  final  time  or  reaching 
a  value  of  infectives  below  an  assigned  minimum,  is  satisfied. 

Solution  of  the  optimum  control  problem  is  obtained  by  an 
iterative  procedure.  First  the  values  of  controls  are  assumed  zero 
and  the  final  time,  which  is  the  time  taken  by  the  disease  to  subside 
without  any  control  being  applied,  is  calculated.  Using  this  time 
as  the  terminal  time  the  system  iterates  to  satisfy  the  cost  immuniz¬ 
ation  criterion,  updating  the  control  vector  at  the  end  of  each 
iteration.  The  computation  stops  when  no  significant  further  improve¬ 
ment  in  the  cost  is  possible;  thus  generating  the  best  control 
strategy  at  the  end. 

The  computer  program  consists  of  a  supervisory  main  program 
and  nine  modular  subroutines.  A  brief  description  of  the  program 
modules  is  as  follows: 

1.  MAIN  PROGRAM:  It  is  the  supervisory  program  which  calls  various 
subroutines,  when  needed,  during  the  computation  of  the  solution. 
Input  to  the  program  is  the  parameters  of  the  disease,  delay  in 
control  (if  anv) ,  solution  termination  criterion  and  penaly  constant. 
Final  control  vectors  are  printed  at  the  end  of  the  execution. 

The  intermediate  results  are  printed  by  the  various  subroutines. 


' 

’ 


Ill 


2.  SUBROUTINE  WGHMLT :  This  subroutine  calculates  the  normally 
distributed  weighting  multipliers  from  the  latent  and  incubation 
periods  of  the  disease.  The  resulting  mul ti pi yi ng-constants  are 
stored  as  two  vectors:  WLP  and  WIP,  at  intervals  of  one  day.  This 
subroutine  is  called  only  once,  i.e.  in  the  beginning. 

3.  SUBROUTINE  RKST:  This  subroutine  solves  the  state  equations, 
using  Runga  Kutta  procedure  for  the  solution  of  differential 
equations.  It  uses  subroutine  RKGS  from  the  SSP  library  and  is 
called  by  both  the  main  program  and  the  search  subroutine  ASRCH. 

This  subroutine  also  evaluates  the  cost  function  J. 

4.  SUBROUTINE  FCT :  This  subroutine  defines  the  state  equations  of  the 
model  for  subroutine  RKGS. 

5.  SUBROUTINE  OUTP:  This  subroutine  stores  and  prints  the  output  of 
subroutine  RKST.  Each  time  RKST  is  called,  the  stored  values  of 
states  are  updated  by  this  subroutine.  Printing  of  the  results  is, 
however,  skipped  if  variable  KEY  is  greater  than  10.  This  is  done 
when  subroutine  RKST  is  called  by  the  search  subroutine  ASRCH. 

6.  SUBROUTINE  RKCST :  This  subroutine  solves  the  costate  equations 
using  the  Runga  Kutta  method.  The  solution  is  carried  out  backwards 
in  time  usinq  subroutine  RKGS.  It  also  evaluates  a  new  gradient 
vector  each  time  it  is  called. 
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7.  SUBROUTINE  FCTCS:  This  subroutine  supplies  the  costate  equations 
to  the  RKGS  subroutine  in  subroutine  RKCST. 

8.  SUBROUTINE  OUTPCS:  This  subroutine  stores  and  prints  the  results 
of  the  solution  of  the  costate  equations. 

9.  SUBROUTINE  ASRCH:  This  subroutine  conducts  a  search  for  the 
multiplying  vector  a  (for  negative  gradient  or  acceleration  steps) 
to  minimize  the  cost  function  J.  The  search  is  conducted  both  for 
a-j  and  in  turn,  using  a  combination  of  bi sectional  search 
procedure  and  interpolation.  It  calls  subroutines  RKST  (with 
printing  of  results  switched  off)  and  interpolation  subroutine 
QDINPL. 

in.  SUBROUTINE  ODINPL:  This  subroutine  interpolates  the  value  of  a 
for  minimum  cost  and  is  frequently  called  in  subroutine  ASRCH. 

A  list  of  the  program  described  above  is  given  in  the 


following  pages. 
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c 

c 

C  **  OPTIMUM  CONTROL  GF  EPIOEMICS  USING  PARALLEL  TANGENT  TECHNIQUES  ** 

C  ******  N.K. GUPTA.  ELECTRICAL  ENGG.  DEPARTMENT  ***** 

C 

C  **  SST(25.200)  IS  THE  STORAGE  APR  AT . STORAGE  IN  ROWS  1-4  IS  FOR  STATES, 
C  5-6  FOR  INFECTIVE  AND  REPORT  RATES,  /-  1  0  FOR  CO- S T A T 1 5  .  I  1  - l 6  FOR 

C  CCNT SOLS .  l 7-20  ARE  FOR  GRADIENT  VECTORS.  10-20  FOR  GRADIENT  FOR 

C  SEARCH, 21-24  ARE  FOR  RIGHT  HAND  SIDE  OF  STATE  EQUATIONS  ANO  25  HAS 

C  STORED  HAMILTONIAN  ** 

C 

C  **  WLP  AND  W[o  ARE  NORMALLY  DISTRIBUTED  WEIGHTAGE  VECTORS  ** 

C 

IMPLICIT  REAL*R( A-H.O-2) 

EXTERNAL  FCt,CUTP,FCTCS,OUTdcS 

D  I  MENS  ION  SST(  25.  20  0)  .  US«CH  (  2  *  20  0  )  ,  WLP  (501  ,’WIP(50>.CST(50).U(2). 
22(200  )  .  BZ(2>  .Y(  20  0  )  ,.DE»Y(5  >  ,PRMT(  5)  ,  AUX<  8.5) 

DOUBLE  PRECISION  LMD3.LMD4 , ACT V,R«VD ,flLMD. BET  A , AK, BK.CK , OK, CRC, PK, 
2PW1,PW2,PW3.PW4,DW5,J0ELAY.TMEND 
INTEGER  ITRN.INT,  INTSTR,  INTCS , MI d, MLP, K EY,  I OEL A Y.  I  SN 
COMMON  SST. USRCH ,WLP,WIP,CST.U.Z,BZ 

COMMON  LMD3,LMD4,ACTV.RMV0,RLMD.BETA,AK.gK,CK.DK.CRC.PK,PWl,PW2. 
2PW3  »PW4 .PW5.U0ELAY .TMEND.  ITRN.INT,  I  NT STR.INTCS.MIP.MLP. KEY,  IDELAY 
901  FO  RM AT (9F10.0) 

905  FORMAT(BIIO) 

911  FORMAT (  1HI//.40X.*  *****  OPTIMUM  SOLUTION  OF  TIME  DELAY  EPIDEMIC  MD 
20EL  ****** > 

912  FORMAT ( I  HO . 1 4X , • THE  DISEASE  HAS  *  4C0NT  ACT  RATE**  *.F5.2.*  INITIA 
2L  SUSCEPTIBLE  POPULATION  0.50  ANO  DISEASE  F  ACT  CR  *.F5.2) 

914  FORMAT!  1MO ,42X ,•  ***  ACTIVE  CONTROL  BECOMES  EFFECTIVE  AFTER*. 14.* 
20AYS  ***• ) 

915  FORMAT (  IH  . 1  OX ,  •  E XPECTED  MEAN  *  LATENT  *  AMD  *INCUBATION*  “ERIODS 
2  ARE  * . 2F6 . 2  »  *  AND  RESPECTIVE  STANDARD  DEVIATIONS*  »2F6.2//> 

916  FORMAT (  l H- ,60X .« COST  CONSTANTS  AR E* // .  4 5X  .  •  CCS T  PER  REPORTED  CASE 

2  =• .F9.2/.45X, 'PER  UNIT  COST  OF  ACTIVE  CONTROL  = 

3* . F9 .2/. 45X , • SQUARE  LAW  COST  OF  ACTIVE  CONTROL  - • ,F9. 2/, 45X • 

4  *  PER  UNIT  COST  OF  PASSIVE  CONTROL  = •  , F9 . 2 / . 4 SX * • S QUA  RE  LAW  CDS 

5T  OF  PASSIVE  CONTROL  = • , F 9 . 2/ , 4 5X .  'PENALTY  FOR  VIOLATING  CONSTR 
6 AINT  S  =*  . F9 .2// ) 

917  FORMAT! IH-.40X,* ******  NORMALLY  DISTRIBUTED  WEIGHTAGE  MULTIPLIERS 
2******* // ) 

918  FORMAT  {  3  {  l  X  .*  DELAY  IN  D  A  Y  S  *  *  *  *  •  _,  2X  ,  l  0  (  3X  .  !  2  .  6  X  )  /  /  .  1  X  .  •  L  .  P  .  M  UL  T  I  PL  I 

2ERS**  *  .  2X.10F11.6//.I  X»*I  . P . VUL T I  PL I ERS *  * •  ,  2X.1 OF  11.  6 /////). 

3  55X  *******  SOLUTION  BEGINS  ♦****•> 

919  FORMAT  (IX.*  PENALTY  WEIGHTAGE  MULTIPLIERS  PWl  FOR  NEGATIVE  I NFE 

2CTIVES=  •  » F8 • 2 • •  PW 2  FOR  RESIDUAL  INFECTIVES=  *.F6.2//> 

920  FORMAT ( IH1/// ,45X, • ITERATION  NUMBER*,I4.*  SOLUTION  WITHOUT  CONTRO 
2L* //) 

921  FORMAT ( IH1 /// ,45X, • I TERATI ON  NUMBER*. 14,*  SOLUTION  WITH  CONTROL*/ 
2/> 

92  2  FORMAT (  1H-, 1  OX, 'DAY'  ,7X.*  REP. RATE*  ,7X.*  REPORTED*  ,5X. •  INFCT .RATE*  . 

2  5X, • INFECTIVES* ,6X. • SUSCPTBLS* ,5X. * ACTV  CNTRL* , 5 X , • PSS V  CNTRL*/) 
925  FORMAT (  1H-,  10X  .*  ***  FINAL  T I ME= * . F B . 2 .  *  *  * *  » ,  1 0 X , *  * *  *  COST=*. 

2  F 1 0 • 6  «  *  ***  *.I0X.****  STORAGE  POS l T I ONS= * . 1 5 , *  ***•) 

935  FOPMAT ( 1 H- . 10X ,* ITERATIONS  HAVE  EXCEEDED  THE  LIMIT*. 16) 

938  FORM  AT (  1  Hi ///  .55X . • ***  FINAL  ACTIVE  CONTROL  ***•) 

939  FORMAT ( 1 H— /// *55X . • **♦  FINAL  PASSIVE  CONTROL  **★») 

940  FORM  AT (  1 H- ./ ,  ( 2X. 1  OF  1 3.6 )  ) 

945  FORMAT (  1HI///.50X. *  ***  ACCELERATION  STEP  FOR  ACTIVE  CONTROL  ***♦) 

946  FORMAT ( 1 H  ///.50X.****  ACCELERATION  STEP  FDR  PASSIVE  CONTROL  ***•> 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


**  READING  VALUES  OF  LATENT  AND  INFECTIOUS  PERIODS,  THEIR  MEAN  VALUES 
AND  STANCARC  DEV  I  AT  I O NS , CONT AC T  RATE  AND  DISEASE  FACTOR.  *  * 

**  ALSO  READING  COST  CONSTANTS  FOR  ACTIVE  AND  PASSIVE  C ONT ROL S ( AK , 3K , 
CK.DKJ.COST  PER  REPORTEO  CASE  CRC,  AND  PENALTY  CCNSTANT  PK .  ** 

**  GRDSUM  IS  THE  FLOOR  LIMIT  ON  THE  SUM  OF  SQUARE  OF  GRADIENTS.  IEXIT 
IS  THE  NUMBER  OF  ITERATIONS  ABOVE  WHICH  ITERATION  SHOULD  STOP  AND 
l DELAY  IS  THE  TIME  DELAY  IN  THE  EFFECTIVENESS  OF  ACTIVE  VACCINE  ** 

RF AD (S.901)  DML3. DM  IP. SOLP, SD IP. BETA ,F ACK ,U ( t ) ,U  C2 ) 

READ (5, 90 1 )  AK , 6K , CK . DK * CRC . PK . GRDSUM 
R F  AC { 5 » 905 )  IEXIT.  IOELAY 

**  Pkl  AMD  PW2  ARE  PENALTY  MULTIPLIERS  FOR  NEGATIVE  AND  RESIDUAL 
INFECTIVFS  RESPECTIVELY.  PW3  IS  THE  VALUE  OF  INITIAL  INFECTIVES. 

P  W  4  AND  P  >  5  APE  CONDITION  CODES  LSED  IN  THE  PROGRAM.  ** 

RE  AD { 5 . 90 1  )  PW 1 .PW2.PW3.PW4 

**  INITIALISING  THE  VARIABLES  TO  ZERO.  ** 

CO  ST=  0 • 0 
TMEND=0 .0 
DO  10  1=1.50 

CST ( I )  =  0. 0 
10  CONTINUE 

DO  20  1=1,25 

DO  20  J-i .ZCC 

IF (  I .LT. I  I .OR. I «GT . 12 )  GOTO  15 
I J=I-10 

SST ( I. J)=U( I J) 

GOTO  20 
15  SST ( I »  J ) =0. 0 

IF (  I .GT .2  )  GO  TO  20 
USRCH (  I , J )=0 .0 
20  CONTINUE 
I  TRN-=  1 

WR  I TE (6.91 1 ) 

WP  I TE (6.912)  BETA  ,FACK 
WRITE(6.915)  OMLP. DM  I P . SDLP . SD I P 
WPITE(6,916)  CRC, AK  .BK.CK.OK.PK 
WPITE(6.919)  PW1.PW2 

**  EVALUATING  AND  STORING  THE  WEIGHTAGE  MULTIPLIERS  FOR  NORMAL 
DISTRIBUTION  OF  LATENT  AND  INFECTIOUS  PERIODS.  ** 

CALL  WGHMLT ( WLP.WIP.DMLP. DM  I P . SDLP. SOI P . BET  A , FACK , M I P , MLP ) 


WRITE! 6.91 7) 

WPITE(6»910)  (I.  1=1.10).  ( WLP (I).l=l«10),(WlP(  I).I-1  .10). 

2(1.1  =  11. 20).  (  WLP  (  X  )  »  I  —  1  1  »20)«(WIP{  I  )  i  1  =  U»20)  » 

3(  1.1  =  21.30).  ( WLP (I  )  .1=21  .30)  , I W IP { I )  .1=21,30) 
WRITE(6,914)  IDELAY 
DO  25  K= l ,200 
USRCH( 1 .K ) =S5T ( 1 1 , K ) 

USRCH(2,K)=SST ( 12. K ) 

25  CONTINUE 
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WRITE(6,920)  ITRN 
WRITE(6,922) 

KEY=! 

CALL  RXST(  XtPSMT,Y,DEHr,FCT,r,UTP,AUX,COST) 

CST(ITRN)=  COST 

*01  TE  (  6  ,  92  5 )  T MEND, COST. INTSTR 

CALL  RKCST(X,PRMT,Y.DE-.Y,FCTCS,OUTPCS,AUX> 

c 

C  **  ONE  DIMENSIONAL  SEARCH  WITH  STEEPEST  DESCENT  FOR  FIRST  STAGE.  “A 

c 

CALL  ASPCH (X ,PPMT , Y , DERY.FCT.CUTP. AUX, COST ) 

c 

C  **  ITERATION  USING  PARALLEL  TANGENTS  TECHNIQUE  BEGINS  *♦ 

c 

35  CONTINUE 

C  ALt_  RKCST  (  X  «PRMT  , Y . DER Y, FC TC S . OUTPCS  ,  AUX  ) 

A=ezt i ;+ezt?> 

IFC A.LT.GRDSUV)  GOTO  60 
KEY=3 

CALL  ASPCHIX.PRMT.Y, DERY.FCT.CUTP, AUX, COST) 

DO  45  V  = I  , I NTSTR 
SST  (  19»M)=SST{  1  1  ,M)— SSTI  13,  M) 

SST(20.M)=SST< 12,M)-SST{14.M) 

45  CONTINUE 

**  W  F  I  TNG  ACCELERATION  STEP  SEARCH  VALUES  ** 

C 

■  FITE ( 6, 945  ) 

WP1TF(5»940)  ( SST ( l 9, N ) ,N=M IP, INTSTO ) 

ic(  o  «  9<*  6  t 

WFITE(6,940>  (SST( 20. N).N=MIP. INTSTR) 

KEY  =  5 

CALL  ASPCH I X.PRMT, Y.DERY.FC  T.OUTP .AUX, COST) 

C 

c  **  ITERATION  TERMINATION  CRITERION.  ** 

c 

IF(  ITRN.GE.IEXIT)  GOTO  55 

GC  TO  35 
55  CONTINUE 

WRITE (6,935)  ITRN 
60  CONTINUE 
C 

c  **  WRITING  FINAL  CONTROL,  ** 

c 

WRI  TE  (  6,93  8) 

WRITE  16, 9*40)  <SST(il.N),N=MlP. INTSTR) 

WPITEt  6,939) 

WRITE(6,940)  (SST(  12. N),N=M IP,  INTSTR) 

STOP 

END 


TOTAL  MEMORY  REQUIREMENTS  00I9PA  BYTES 
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c 

c 

c 

SUBROUTINE  V/GHMLT  (  WLP  .  W  I  P,  OMLP,  DMI  P,  SOL°,  SOIP.OETA.FACKtMIP.Mi_P) 
C 

C  **  THIS  SUBROUTINE  CALCULATES  THE  NORMALLY  DISTRIBUTED  MULTIPLIERS 
C  FOP  THE  CALCULATION  OF  ACTIVE  AND  REPORTED  CASES.  ** 

C 

IMPLICIT  PEAL*at A-H.O-Z) 

DIMENSION  WLPI50 ) . W  IP{ 50 ) , ST ( 50 ) 

DO  10  1=1,50 

WLP  {  1  ) =0 .00000 
W  IP  ( I ) =0 .00000 
ST( I )  =0. 0 

10  CONTINUE 
C 

C  ***L0ADINC.  THE  MULTIPLICATION  ARR  AYS**** 

C 

VLP=2.0*0MLP+I 
MI P=2»  0*DMI P+1 
DO  45  1=1.2 

IF(I-i)  20.20.22 
20  K=MLP 
SD=  SDL.P 
DME  AN=OMLP 
GO  TO  25 
22  K = M  IP 

SD=SD I P 
DMEAN=OMIP 

25  C=FACK48ETA/(SD*DSQRT(5.2831D0) ) 

DO  45  J= 1 • K 
T  =  J—  1 
C 

C  **  CALCULATING  THE  NORMALLY  DISTRIBUTED  CONSTANT  MULTIPLIERS.  ** 

C 

ST( J)=C*OEXP( (  —  0.5*(D  MEAN— T ) * *2 ) / ( S D* *  2 > ) 

C 

C  **  STORING  THE  MULTIPLIER  CONSTANTS  IN  PROPER  ORDER.  ** 

c 

M=M I P  + 1  —  J 
IF(I-l)  30.30,35 
30  WLPI M) =ST( J) 

GO  TO  45 
35  WIP(M)=ST< J) 

45  CONTINUE 
RETURN 
END 


TOTAL  MEMORY  REQUIREMENTS  00062C  BYTES 
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c 

c 

c 

subroutine  RKST (X  .prmt .  y  .  de  r  y »FC“.OUTP»  au  x « cc s T ; 

♦♦this  subroutine  integrates  the  state  equations  using  runga-kutta 
method  and  subroutine  RNGS  from  ssp  library.  ♦  «, 

IMPLICIT  REAL*S( A-H.O-Z ) 

DIMENSION  SST ( 25.  200) .USRCHI 2.  200 ) «  #L3 (50).WIP(50).CST(50)  .UC2)  . 
2ZC  200)  .BZ( 2)  .Y(?00)  ,0E  *Y{5)  , PRMT (5)  .  AUX (8.5 ) 

DOUBLc.  PR  EC  T  S  I  ON  LMD3.L*0<.  ,aCTV»RMVD«RLM.d»8CTA.4K«3K»CK»OK»  CRC  «  PK  » 
2  P  '*•  1  .  P  W2  iP».J  i  PW4  .PW5»'J  DELAY  •  T  MEND 
INTEGER  I TRN . 1  NT . i NT  STR , I NTCS .Ml P, VLP. KEY . I  DELAY . I SN 

common  sst .usrch.wlp.wtp.cst.u.z.bz 

COMMON  LMD3.LM04 . ACTV . PMYD . RLMD. BET A , AK  .  BK • CK . DK.CRC. PK • Prfi *  PW2. 

2P  W3,  PW4,  PW5.UI)cU  Y  ,  T  '"END,  I  TKN,  I  NT  ,  I  NTSTR  ,  I  NTCS  .MI  P.MLP  ,KEY  ,  I  DELAY 

921  FORMAT  i  IM1///.45X.  »  ITERATION  NUMBER •  . I  4 •  •  SOLUTION  WITH  CONTROL  * / 
2/) 

922  FORM ATI  IH-,  1  OX. • DAY •  •  7X  .  *  PEP. PATE'  .  7X,  •  REPORTED*  ,5X.*  I  NFC T. RATE*  . 

2  5X,  •  INFEC  TI  VES  •  ,  (>X.  •  SUSCPTELS  •  ,  5X.  •  ACTV  CNTRL  *  *  5X,  •  PS5V  CNTRl*/) 

C 

C  ♦*  INITIALISING  PARAMETERS  AND  STATES.  ** 

C 

PPMTI 1 )~-l . 0 
PPMT(2)=0.0 
PRMT (3)=0.25 
PRMT ( 4) = 1 . OD-5 
Y{  1  )=0.5 
YI2)=PW3 
Y(3)=0.0 
Y(4)=0.0 
I  N  T  =  0 

SST( 1 ,MIP)-Y( 1 ) 

SST ( 2.MIP) =Y ( 2  > 

c 

C  ♦ ♦ I NT  REPRESENTS  THE  T I Vf  IN  DAYS  ** 

C  ♦♦  UPDATING  PARAMETERS  FOR  EACH  INTEGRATION  INTERVAL.  *♦ 

c 

10  CONTINUE 
INT=INT*1 

PRMT ( 1 )=PRMT{ 1  )+l  . 

PRMT  C2)=PRMT ( 2) +1  . 

DO  20  J= l .4 
DER Y C  J ) =  0.25 
20  CONTINUE 
ACTV=0.0 
RMVD=0.  0 

c 

C  ♦♦  CALCULATING  ACTIVE  AND  REPORT  RATES  FOR  GIVEN  DAY  *♦ 

C 

DO  30  J—  I  .  M I P 

M= I  NT* J-I 

Y  2— SST (  t ,M ) *SST{ 2 . M) 

ACTV= ACT V4WLPI J)*Y2 

RMVD=RMVC+WIP(J)*Y2 
30  CONTINUE 
C 

C  ♦♦  PRESENT  VALUES  ARE  STORED  AFTER  MlP  PLACES  •* 

C 

N— 1 NT+M  I P- I 
NCCL  AY  =  N—  I  DELAY 
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S3T(5,N>=ACTV 
SST l 6  »N ) =R V VD 
U(  1  )=USRCH<  1  ,N) 

UCEL  AY  =  LSRCI-  (  1  «  NOEL  AY  ) 

U(  2  )  =  USRCt-  (  2  ,N  ) 

CALL  DRK  GS  < PRVT  ,Y  .DERY . 4 • IHLF.FCT .OUT3,  AUX  5 
IF(KEY-l)  35.35.40 

35  IF{Y(2).GE.PW3.A.ND«FRmT{2).LT.12«0DI>  GOTO  10 
C 

C  *♦  T  MEND  IS  THE  FINAL  TJV£,  l NT  END  IS  THE  TOTAL  NUMBER  OF  INTEGRATION 
C  INTERVALS  AND  I  NT  ST  R  IS  THE  NUMBER  OF  ST  OP A  GF  POSITIONS  NEEDED.  ** 

C 

TVEND=RRMT (2 ) +25. 

INTEND^ INT+25 
I  NT S T R—  M I P  + INTEND 
K  E  Y-=  2 

40  IF( INT.NE. INTEND)  GOTO  10 

**  CALCULATING  THE  COST  DEPENDENT  ON  FINAL  VALUES  OF  STATE.  ** 

c 

CCST=CRC*(PW2*Y(2)+Y(3))+Y(4) 

RETURN 

ENO 


TOTAL  MEMORY  REQUIREMENTS  00064E  BYTES 


C 

C 

c 

SUBROUTINE  FCT (X.Y.DERY) 

c 

C  **  THIS  SUBROUTINE  GIVES  STATE  EQUATIONS  WITH  CONTROL  ♦* 

c 

IMPLICIT  REAL*9 < A-H.O-Z ) 

DIMENSION  SST ( 25. 2  00  ,USRCH(2.200)  . WLP ( 50  >  .WJP(50).CST (50  ).U(2)  . 

2  Z { 20  0 )  • BZ ( 2 )  .  Y (  20  0)  •  DERY ( 5 )  . P  R  M  T 15)  . AUX ( 9  *  5 ) 

DOUBLE  PRECISION  LMD3ILMD4. ACTVi RMVO. RL^D. BETA . AK, BK.CK «0< . CSC ,PK. 

2PW1  .PW2.PVv3.PW4  .PW5.UDELAY,  TMEND 
INTEGER  ITRN,  I  NT.  I NT  ST R.  I NTCS • m I P >  MLR «  KEY .  IOELAY.ISN 
COMMON  SST , USRCH .WLP.WIP.CST.U.Z.BZ 

COMMON  LMD  3  »  L  WD4 , ACT  V  »  R  MVO  *  RLMD  »  BET  A  »  AK  «BK  »  CK . DK  »  CRC  »  PK . PW 1  . PW  2 . 

2P W 3.P W 4.P W 5. U DELAY. T MEND.  ITRN.  INT .  InTSTTi  INTCS.MIP.MLP.KEY.  I DE  LAY 
C 

C  **  HY2  IS  THE  CIRAC  DELTA  FUNCTION  FOP  NEGATIVE  Y<2>  PENALTY  MULTIPLIER  ** 

c 

HY2=0. 0 

DEPY  (  1  )•=  -6FTA*Y(1)*Y(2  )  — UDELAY  4  Y  (  l)-U( 2) *Y< 1 ) 

OERY (  2  )  =  ACTV-PMVO-U< 2 ) *Y( 2 ) 

DERY(3>=  RMVD 

IF( Y( 2 ) .LT .0 .0  )  HY2= 1.0 

DERY (4)=PK+PWl WHY 2*Y(2)**2+( AK  +  SK *U  (  l)  )*UC  1  >+(CK  +DK  *U(  2  )  )*U(2) 

RETURN 

END 


TOTAL  MEMORY  REQUIREMENTS  000374  BYTES 
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c 

c 

c 


c 

c 

c 

c 


**  this  subroutine  stores  the  current 
the  results.  ** 


VALUES  OF  STATE  AND  PRINTS 


IMPLICIT  REAL*B ( A— H, O— 2 ) 


c 

c 

c 

c 

c 

c 

c 


910  FORMAT (  1H  . F 1 3 • 2 , BF I  5 • ft ) 

I  F  (  X.  NE.  0 . 0  .  AND  .  X  .  NE.  PRMT  (2  )  )  GO  TO  55 
IF(X.EQ.O.O)  GOTO  40 
I  K= I  NT  +M  I  P 
DO  1  0  J-l  ,  4 

10  SST( J, IK )=  Y( J ) 

IF( KEY. GE. 1 0 )  GO  TO  55 

DO  20  J= 1 , 4 

JN=20+J 

20  SST ( JN, IK > =DERY ( J ) ' 

IF( Y( 2) .LT . 0.0 )  GO  TO  50 
IF< INT/5*5.NE. INT )  GO  TO  55 
GOTO  50 

40  IF ( KEY .GE . i 0 )  GO  TO  55 
DO  45  J= 1 , 4 
JN= J+20 

45  SST( JN.MIP)=DERY( J) 

50  WRITEI6.910)  X , RMVD, Y ( 3 ) . DER Y { 2 ) . Y I  2 )  .  Y  (  l ) . U(  1  ) . U (  2 ) 

55  CONTINUE 
RETURN 
END 


TOTAL  MEMORY  REQUIREMENTS  0C0544  BYTES 


SUBROUT  I NE  RKCST ( X . PRMT. Y. DERY . FCT CS * OU TP CS . AUX ) 

♦♦THIS  SUBROUTINE  INTEGTRATES  BACKWARDS  THE  COSTATE  EQUATIONS  USING 
SUBROUTINE  R KGS  FROM  THE  SSP  LIBRARY.  ** 

IMPLICIT  REAL*B< A-H.O-Z) 

DIMENSION  SST { 25 .200)  .USRCH (2.200).WLP(50).WIP(50).CST(50).U(2). 
2Z(200)»BZ( 2)  «Y(  200)  .DERYl 5)  .PRMT (5)  .AUX(8.5) 

D  OUBLE  PRECISION  LMD3.LM04. ACTV.RMVD.RLMD.BETA .AK.3K.CK.DK.CRC.PK. 
2P W1.PW2.PW3. PW4 ,PW5. UDELAY. TMEND 

I  NT EGER  ITRN. INT.  I NTSTR . INTCS.MIP.MLP.KEY. IDELAY . ISN 
CCMMON  SST.USRCH.WLP.W IP.CST.U.Z.BZ 

COMMON  LMD3.LMD4.ACTV.RMVD.RLMD.3ETA.AK.BK.CK.DK.CRC.PK.PWl.PW2. 

2  pw  3  .  py«  4  «PW5, UDELAY. TMEND.  ITPN,  INT  .INTSTR.INTCS.MIP.mlR.KEY.  IDELAY 


c 

c 

c 

c 

c 

c 

c 
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915  FORMAT (  1H1 /// .24X  .  •  ***  SOLU  T I  CM  OF  COSTATE  EQUATIONS  AND  CALCULATE 
20  VALUES  OF  HAMILTONIAN  AND  GRADIENT  VECTORS  ***•) 

920  FORMAT  (  1H— .  10X  ,  ‘DAY  *  «  8  X  «  ' LEMOA-1  •  «  8  X  »  'LEM DA- 2*  . 8 X .  • LE MO A- 3 •  , OX . 

2  *  LF  M0A-4 •  »  4X  .  •  HAMILTONIAN'  ,5X,  • GR AO  I  ENT- 1  •  ,  5X  , 'GRADIENT  — 2*/) 

950  FORMAT (  1H0 «  10X ***  GRADIENT  I NTEGER AL S* *  *  '.4F15.6) 

C 

C  **  INITIALISING  The  PARAMETERS  FOR  BACKWARDS  INTEGRATION  OF  COSTATES. 
C  INTCS  IS  THE  STORAGE  INDEX  STARTING  BACKWAROS  FROM  LAST  STORED 

C  VALUE.  +* 

c 

INTCS= INTSTR+ 1 

I  NT  =  0 

Y(  l  )=0.0 

Y ( 2 )=PW2*CRC 

LMD3=CRC 

LMD4=  1  . 

PRMT (  1  )=TMEND+ 1 .0 
PR  MT ( 2 l=TMEND 
PRM  T ( 3 ) — — 0.25 
PRMT  I  4  )  =  1  •  0  0-5 
W  R I T€ (6*915) 

WRITE(6.920) 

SST (7, INTSTR  >=Y (  1  ) 

SST( 8. I NTSTR )=Y (2 ) 

C 

C  **  UPDATING  THE  PARAMETERS  FOR  BACKWARD  INTEGRATION  INTERVALS.  ** 

C 

10  CONTINUE 

I  N  T  C  S=  I  NT  C  S—  1 
PRMT II) =PRMT ( l )-l .0 
PPMT(2)=PRMT(2>-1.0 
DERYI 1 >=0.5 
DERYI 2)=0.5 
RLMD=0 .0 

IFIINT.EQ.O)  GOTO  20 
IF ( I NT.GT. MI P)  INT=MIP 
DO  30  1= 1  * INT 

M  W=  M  I  P—  I  +  1 
MCS  = I  NT CS ♦ 1-1 

RLMD  =  PLMD+ WLPI  MW)*SST(8.MCS)-WIP(MW)*(SST{8.  MCS  )  -LMD3  ) 

30  CONTINUE 
20  I NT= I  NT  +  1 

CALL  DR KGS (PRMT .Y.DERY.2. I HLF . FCTCS . OUT PCS , AUX) 

I F ( PRMT (2J .GT .0. )  GO  TO  10 

DO  40  L=MIP. INTSTR 
SST(  19  »L )=  —SST (  1 7 . L ) 

SST<20.L>=  —SST ( 1 8 . L ) 

40  CONTINUE 
C 

C  **  CALCULATING  INTEGRAL  OF  SQUARE  OF  GRADIENTS  ** 

c 

DO  50  1=1*200 
Y  <  I  )  =  0.0 
50  2<I)=0.0 

K=I NTSTR+l-MIP 
DO  70  L= 1*2 
J=L+18 
DO  60  1=1 .K 
ML= I  —  1 ♦MIP 
Y {  I ) =  SST (  J  *  ML ) **2 
60  CONTINUE 

CALL  DOSF (  1 . 00* Y . Z. K  ) 

BZ(L)=Z(K> 

70  CONTINUE 

WPITEI6.950)  BZflJ.B Z(2) 

RETURN 

ENO 

TOTAL  Mf MfRY  REQUIREMENTS  000700  BYTES 
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c 

c 

c 

SUBROUTINE  FCTCSI X, Y , OERY ) 

c 

c  **  this  subroutine  gives  costate  equations  ** 

c 

implicit  real*b{ a-h.o-z ) 

DIMENSION  SST<  25,2  00)  .USRCH (2,200). *LP<50> , W I P( 50 ) « CST { 50 ) , U ( 2  > , 

2Z (200 ) • 92 ( 2)  ,Y<  20  0 ) ,DERY( 5) .PRMTI5)  . AUX( 8,5) 

DOUBLE  PRECI  SION  LMD3 ,LMD4.  ACTV,  PMVO,  RL MD. BET A . A K . BK , CK  ,DK ,CRC ,PK. 
2PW1  ,P»2,PVk2  ,PW4  ,PW5,UDELAY,  TMENO 
INTEGER  IT  RN. I  NT ,  I  NT  STR ,  I N T C S , M I P , MLR , KEY ,  I  DELAY , I SN 
CCFMON  SST .USRCH , WLP.W I P, CST, U. Z ,BZ 

COMMON  LMD3,LMD4,  ACTV. PM VP,P.LMD, BETA, AK  ,BK,CK.DK,CRC.PK,Prfl  ,PW2. 
2PW3.PW4 . PW5.U0ELAY.TMEND.  ITRN,  I  NT, I  NT STR,  INTCS.MIP. MLP .KEY , IDELAY 
HY2=0.  0 

MDL  AY  = I  NT CS- IDELAY 

DERYC  1  )=  Y(  1  )*(SST(  1  1 . MDL  Av  )+BETA*SST(2.INTCS)+SSTtl2.INTCS))- 
2  RLMD* SST ( 2 . I NTCS) 

IF(SST( 2, INTCS)  ,LT  .0.0 )  HY2=I.O 

DER Y ( 2 ) =  Y(1)*8ETA*SST(1,I NTCS ) -2 . 0 *LMD4*PK *SST ( 2, INTCS ) *HY2*PWl 
2  —  RLMD  ASST  (  l . I NTCS) *SST(  12,INTCS)*Y(2> 

RETURN 

END 

TOTAL  MEMORY  REQUIREMENTS  000442  BYTES 


C 

c 

c 

SUBROUTINE  OUTPCS(X«Y»DERY,IHLF,NDIM«PRMT) 

C 

c  **  THIS  SUBROUTINE  STORES  THE  COSTATE  VALUES  ** 

C 

IMPLICIT  RE AL*0 ( A-H, 0-Z ) 

DIMENSION  SST(25.200) ,USRCH(2»200)»WLP(50)»WIP(50)»CST(50)»U(2)« 
2ZI200)  ,9Z(2)  ,  Y (  200)  ,0EPY(5)  ,PRMTI5)  »AUXt8»5) 

DOUBLE  PRECISION  LMD3.LM04, ACTV ,RMVD, RLMD, BETA, AK,BK,CK.D<»  CRC , PK , 
2PW  1  ,P  *2  »  P*3  ,  PW4  ,P#5  ,  UDELAY,  TMEND 

I  NT EGER  ITRN.INT.INTSTR,  INTCS.  MIP,MLP , KEY , IDELAY,  ISN 
COMMON  SST , USR  CH , WLP, W IP, CST ,U,Z.BZ 

COMMON  LMD3,LMD4»ACTV,RMVD, RLMD, BETA, AK.BK.CK.DK, CRC. PK.PW1.PW2, 

2P1*3.PW4,PW5»UDELAY  ,  TMEND  .  I  TRN,  I  NT*  I  NT  STR.  I  NTCS  .  M I  P  ,  MLP,  KE  Y  ,  IDELAY 

910  FORMATt  1H  ,F  1  3. 2 . BF 1 5. 6 ) 

IF( X ,NE . TMEND. AND . X. NE . PRMT ( 2 ) )  GOTO  55 

IF! X.EQ. TMEND )  GOTO  20 
C 

c  **  STORING  the  VALUES  OF  COSTATE  ** 

C 

L- INTCS- 1 

SST ( 7.L )=Y ( 1 ) 

SST(8,L)=Y(2) 

GOTO  30 
20  L= 1NT  CS 
30  M=L*IDELAY 


n  n  n 


-1  22- 


**  CALCULATING  THE  HAMILTONIAN  ** 

SST( 25 ,L) =SST( 21  .L) *SST ( 7 ,L ) +SST ( 22. L>  *SST ( e.L  >+SST(23. L ) *LMD3  ♦ 
2  SST! 24 ,L ) *LMD4 

c 

c  **  CALCULATING  AND  STORING  NEW  GRADIENT  ***** 

c 

SSTC 17,L)=-SST ( 7 • M) *SST| l.M)*LMO4*(AK+2.0*BK*SST( 1 1 ,L  ) ) 

SST  {  I  a  .L  )  =  LMD4  *  (  CK  +  2 .  0*l>K*SST  (12. L)  )-SST!8,L)*SST(  2.L) 

2  -SST(7,L)*SST( l.L) 

IF(X.EQ.TMENO)  GOTO  50 
IE ( INTCS/5*5.NE . I NTCS)  GOTO  55 

50  WRITE(6,910)  X,Y(l).Y(2).LMD3,LM04,SST(25.L).SST(17,L),SST(18.L) 
55  CONTINUE 
RETURN 
END 


TOTAL  MEMORY  REQUIREMENTS  0006FC  BYTES 


C 

c 

c 

c 

c 

c 


SU8R0UT  INE  A  SRCH ( X . PRM T , Y . DER Y . FC  T, 0 UT P . AUX , CO  ST ) 

#*  THIS  SUBROUTINE  PERFORMS  ONE  DIMENSIONAL  SEARCH  ON  ALPHA  ** 
IMPLICIT  RE AL*8 ( A-H, 0-2 > 

DIMENSION  SST (25. 2 00) « USRCH ! 2 » 20 0 ) .WLP(50).WIP(50).CST(50),UI2) . 

2  2(200) »RZ( 2) »Y(  200) .OERYI 5) • PRM  T ( 5)  . AUX I  8 .5 ) 

DOUBLE  PRECISION  LMD3.LMD4, ACTV.RMVD.RLMD. BETA .AK.BK.CK  »DK .CRC .P<. 
2PW1  ,PW2  . PW3  .PW4.PW5  .  UOELA  Y  •  TMEND 

INTEGER  ITRN. INT.INTSTR. INTCS.MIP.MLP.KEY.I DELAY . I SN 
COMMON  SST • USRCH . WLP. W I P.CST .U.Z.BZ 

COMMON  LMD3.LMD4.ACTV.RMVD.RLMO.BETA.AK.8K.CK.DK.CRC.PK.P'*'  1»DW2. 

2PW3.PW4.PW5. UDELA  Y . T  ME  NO  «  I TRN. I  NT . I NTSTR . INTCS .MIP. MLP .KEY,  IDELAY 
)15  FORMAT!  1H0.10X.  *  *  *  *  CALCULATED  VALUE  OF  ALPHA=  *, F 1 0 • 6 .****•« 1 0 X , 

2  •***  MINIMUM  COST=  •  »F10.6»*  ***•) 

)16  FORMAT ( 1H1  »33X»  *  ****  ONE  DIMENSIONAL  SEARCH  ON  CONTROL  U*  .  I  1  . 

2  •  BEGINS  ****•//»  4  OX . 1  SEARCH  NO. ’.6X« • ALPHA '  »11X» ’COST*/) 

320  FORM  AT (  1 H  *40X»I5. 2F 15.6) 

?2l  FORMAT ( 1 H  ///.45X,  •  ITERATION  NUMBER*.  14.*  SOLUTION  WITH  CONTROL* 

322  FORMAT! 1H-.1 OX.* DAY* , 7 X . • RE P . R ATE • ,7X ,* REPORTED* .5X.* INFCT.RATE* . 

2  5X .* INFECT  IVES*  . 6X .* SUSCPTBLS* ,5X,  ■ ACT  V  CN TRL *  .  5 X , • P SS V  C  N  T  RL  * / ) 

323  FORMAT! 1H  .40X.*ONLY  ACTIVE  CONTROL  HAS  BEEN  SEARCHED  FOR  MINIMUM 
2COST • ) 

324  FORMAT! 1H  ,40X.*ONLY  PASSIVE  CONTROL  HAS  BEEN  SEARCHED  FOR  MINIMUM 
2  COST* ) 

325  FORMAT (  1 H— , l OX , ' ***  FINAL  T I ME=  *  ,F8.2.*  * * * • . 1 0 X . • * *♦  COST- *  . 

_  FIO  6  •  ***  *  .  10X  .  *  ***  STORAGE  PO S I T I O NS= * » I  5 ,  •  ***•> 

FORMAT  C1H  .«2X..8CTM  CONTROLS  HAVE  BEEN  SEARCHED  FOR  MINIMUM  COST- 

2) 

330  FCRMAT11H  , 26X • •  I  NT ERPOL ATEO  ALPHA  ) 

ICHK=0 
ITRN=  ITRN  + 1 
I  SN=  1 

5  CONTINUE 
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c 

c 

c 


c 

c 

c 


Ke  Y  =  KEY  ♦  1  o 
I  CM K=  I  C HK4- 1 
ALPHA -0. 0 
PW4=  1  .0 

**  I SRCH  IS  SEARCH  NUMBER  AND  I  TEST  TESTS  VALUES  FOR  INTERPOLATION  ** 

I SRCH=0 
I TEST=0 

WRITE! 6,91 6)  ISN 

WRITE(6,920)  ISRCH, ALPHA, COST 

K 1=1 SN 

K  2- I SN  +  10 

K3= I SN+ 1 2 

K4=I SN+  1  4 

K5=ISNM8 

I  FLAG=0 

DALPHA=0.0016 

CST 1 =0 , 0 

CST2=0 .0 

CST  3=0 , 0 

A  LP 1  =  0  • 0 

ALP2=0 • 0 

ALP 3=0  «  0 

10  I  TE  ST= I TEST+1 
CSTM IN=COST 
ALPMIN= ALPHA 
C  ST 1 =CST 2 
CST  2=CST  3 
CST3=COST 
ALP  1=ALP2 
ALP2=  AL  °3 
alp3= alpha 
15  ISRCH=I SRCH+1 
PW5=1 .0 

IF{ ISRCH. GT. 10 )  GOTO  35 
IB  IF( ISRCH. GT.l. AND. IFLAG.EQ.O)  D AL PH A= ALPH A # P W5 

ALPHA=ALPM IN+OALPHA 

+  *  UPDATING  THE  CONTROL  WITH  CALCULATED  AL  °HA  FOR  NEXT  TRIAL.  ** 

DO  20  K=MIP,INTSTR 

U SRCH I K 1 ,K )=  SST ( K2 , K) +ALPHA*SST (K5.KI 
IF  t USRCH ( K 1 , K ) • LT .0.0)  USRCH(K 1  ,K )=0 . 0 

20  CONTINUE 

CALL  RK  ST ( X.PRMT , Y, DER Y.FCT, OUTP , AU X , CO  ST } 

WPITE(6,920)  ISRCH, ALPHA, COST 
IF(COST.LT.CSTMIN)  GOTO  10 
A ALP=DABS ( ALPHA ) 

IF! AALP.LT . 1  .0-5 )  GOTO  28 

IF ( 1 FLAG.EO. 0. AND . T TEST. GT .3 )  GOTO  25 

IF(  ITEST . G E • 2 )  GOTO  30 

D ALPH A=D ALPH A/4 .0 

I FL  A  G—  l  OIT 

GOTO  15 

25  D ALPH A  =ALP3— ALP2 
I FL  AG= 1 00 
GOTO  15 

28  IFIPW4.LT. 0. DO)  GOTO  45 
PW4=— PW  4 
PW5=-PW5 
I FL  AG— 0 
I  S  R  C  H=  1  S  R  C  H  ♦  1 
GOTO  18 


u  u 
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**  INTERPOLATION  RUN  FOR  ALPHA  ** 

30  C  ST  1 =C  ST  2 
C  ST  2=C  S  T  3 
CST  3=C0ST 
ALP  I = ALP2 
A  LP  2=  AL  P  3 
ALP3= ALPHA 

I SRCH= ISRCH+ 1 

35  CALL  GDI NPL ( ALP1 . ALP2 . ALP3, CST1 » CST2, CST3, ALPHA ) 
WRI TE<  6,930  ) 

DO  40  K=MIP,INTSTR 

U  SRCH  (  K1.K)=  SST  (  K2  »  K  )  ♦ ALPH A*SST ( K5 , K ) 

IF ( USRCH ( K1,K).LT.O.O)  USRCH(KI,K>=0.0 
40  CONTINUE 

CALL  RKSTt  X.PRMT, Y . D ER Y  ,  FC T , OU TP , AU X , COST ) 
WRITE(6,920)  I  SR C H , AL PHA , CO S T 
45  I F ( CSTM IN. LT .COST )  ALPHA=ALPMIN 

C 

c  **  STORING  NEW  CONTROL  FOR  NEXT  ITERATION.  ** 

C 

DO  70  K=MIP,INTSTR 
I F ( KEY— 13)  55,55,65 

55  SST(K3,K)=SST(K4.K) 

SST(K4«K)=SST(K2«K) 

65  SSTIK2.K )  =  SST { K2  »  K  ) ALPHA  *SST ( K5  »K ) 

IF(SST(K2,K)  .LT.O.O)  SSTCK2,K)=0.0 
USRCHtKl ,K )=SST(K2,K> 

70  CONTINUE 

KFY=KEY-IO 

WR I TE (6,921)  ITRN 

|F(ICHK.E0»1  ,AND.  ISN.EQ.l  )  WRITfct 6,923) 

I  F (  ICHK.EO.l  .AND.  ISN.EQ.2)  WRITE(6«924) 

IF(  ICH< .EO  .2  )  WRITE! 6, 926 ) 

WRITE(6,922) 

CALL  RKSTt  X.PRMT»Y,DERY*FCT»OUTP»AUX«COST) 

C  ST ( I TRN ) =  COST 

WRITE(6,925)  T MEND . COST , I NT STR 
I  ST=1 

I F (  ISN.EQ.2)  IST=- 1 

I SN=I SN* 1ST 

I  Ft ICHK .EQ . 1  >  GOTO  5 

I SN= I SN— 1ST 

RETURN 

END 


TOTAL  MEMORY  REQUIREMENTS  OOOE44  BYTES 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  QDINPLt X1.X2.X3,Y1.Y2.Y3,Z> 


**  this  subroutine  interpolates  values  of  alpha 

SERIAL  VALUES  OF  ALPHA  AND  COSTS.  XI, X2  AND 
ARRANGED  IN  THE  ASCENDING  ORDER  OF  ALPHA  AND 
CORRESPOND  I NG  COSTS,  Y3  3EING  THE  FIRST  VALUE 


MULT  I DL I F R 
X  3  ARE  THE 
Yl.  Y 2,  Y  3 
when  COST 


FROM  THREE 
ALPHAS 
ARE  THE 
STARTS 


INCREASING,  Z  IS  THE  INTERPOLATED  VALUE  OF  ALraH  A .  ** 

ORDER  OF  ALPHA  AND  Y l . Y2  ANO  Y3  ARE  CORRESPONDING  COSTS  Y3  BEING  FIRST 
VALUE  ft  HE  N  COST  STARTS  INCREASING.  Z  IS  INTERPOLATED  VALUE  OF  ALPHA  ** 


IMPLICIT  REAL*8( A-H.O-Z) 

910  FORMAT ( 45X ,• ALPHA  VALUES  DO  NOT  SATISFY  INTERPOLATION  CONDITION*) 
915  FORMAT ( 45X ,* DFNOMI NATOR  =  0.0,  INTERPOLATION  NOT  EXECUTED*) 
X23=X2-X3 
X31=X3— X 1 
X 1 2=X 1  —  X2 
XS23=X 2**2-X3**2 
XS31=X3**2-X 1**2 
XS  1  2=X 1 **2  — X2**2 
XD=2.0  * ( X23*Y l +X3I *Y2  +  XI  2  *Y3  ) 

AXD=DABS( XD) 

IFIAXD.LT .  1 .00-7)  GOTO  15 
D=XO/(X23*X31*X12 ) 

IF(O.GT.O.O)  GO  TO  10 
XN=XS23*Yl+XS31*Y2f XS12*Y3 
Z=XN/XO 
GO  TO  20 
10  WPITEI6.910) 

GOTO  20 
IS  I  TE  (5*915) 

Z  =  X2 

20  CONTINUE 
RETURN 
END 


TOTAL  MEMORY  REQUIREMENTS  0003E0  BYTES 
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